Initial commit from Polytechnique Montreal

1 files changed, 380 insertions, 0 deletions

diff --git a/Dragon/src/EVOSAT.f b/Dragon/src/EVOSAT.f
new file mode 100644
index 0000000..b650be8
--- /dev/null
+++ b/Dragon/src/EVOSAT.f
@@ -0,0 +1,380 @@
+*DECK EVOSAT
+      SUBROUTINE EVOSAT(IMPX,MAXA,MAXB,MAXY,LOGY,NSAT,NVAR,KSAT,YST1,
+     1 YSAT,MU1,IMA,NSUPF,NFISS,IDIRAC,KFISS,YSF,ADPL,BDPL,NSUPFG)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Lumping of the depletion matrix, fission yields, sources and initial
+* conditions to take into account the saturation of depleting nuclides.
+*
+*Copyright:
+* Copyright (C) 2002 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
+* IMPX    print parameter.
+* MAXA    first dimension of matrices ADPL and AGAR.
+* MAXB    first dimension of matrices BDPL, IMA and MU1.
+* MAXY    second dimension of matrix YSF.
+* LOGY    number of passes through EVOSAT:
+*         first pass updates YSAT and YST1;
+*         second pass does not update YSAT and YST1.
+* NSAT    number of saturating nuclides.
+* NVAR    number of nuclides in the complete depletion chain.
+* KSAT    position in chain of the saturating nuclides.
+* NFISS   number of fissile isotopes producing fission products.
+* IDIRAC  saturation model flag (=1 to use Dirac function contributions
+*         in the saturating nuclide number densities).
+* MU1     position of each diagonal element in vector ADPL.
+* IMA     position of the first non-zero column element in vector ADPL.
+* NSUPF   number of depleting fission products.
+* KFISS   position in chain of the fissile isotopes.
+* YSF     product of the fission yields and fission rates.
+* ADPL    depletion matrix.
+* BDPL    depletion source.
+*
+*Parameters: input/output
+* YST1    number densities for all isotopes as input and of
+*         the non-saturated isotopes as output.
+*
+*Parameters: output
+* NSUPFG  number of lumped depleting fission products.
+* YSAT    number densities of the saturating isotopes.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IMPX,MAXA,MAXB,MAXY,LOGY,NSAT,NVAR,KSAT(NSAT),MU1(MAXB),
+     1 IMA(MAXB),NSUPF,NFISS,IDIRAC,KFISS(NFISS),NSUPFG
+      REAL YST1(NVAR),YSAT(NSAT),YSF(NFISS,MAXY,LOGY),ADPL(MAXA,LOGY),
+     1 BDPL(MAXB,LOGY)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(EPS=1.0E-5)
+      CHARACTER HSMG*131
+      LOGICAL LTEST
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: KEV,MGAR,IGAR
+      REAL, ALLOCATABLE, DIMENSION(:) :: YSTG,AGAR,BGAR,GAR
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: A22,YSFG
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: A21,A12
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(KEV(NVAR),MGAR(NVAR-NSAT),IGAR(NVAR-NSAT))
+      ALLOCATE(YSTG(NVAR-NSAT),A22(NSAT,NSAT),A21(NSAT,NVAR-NSAT,LOGY),
+     1 A12(NVAR-NSAT,NSAT,LOGY),AGAR(MAXA),BGAR(NVAR-NSAT),
+     2 YSFG(NFISS,NSUPF),GAR(NSAT))
+*
+      NSUPL=NVAR-NSUPF
+      I0=0
+      DO 40 I=1,NVAR
+      DO 10 II=1,NSAT
+      IF(I.EQ.KSAT(II)) GO TO 20
+   10 CONTINUE
+      I0=I0+1
+      KEV(I)=I0
+      GO TO 40
+   20 DO 25 L=1,LOGY
+      IF(ADPL(MU1(I),L).EQ.0.0) CALL XABORT('EVOSAT: ZERO DIAGONAL COM'
+     1 //'PONENT FOR A SATURATING ISOTOPE.')
+   25 CONTINUE
+      DO 30 II=1,NFISS
+      IF(I.EQ.KFISS(II)) CALL XABORT('EVOSAT: A FISSILE ISOTOPE IS SAT'
+     1 //'URATING.')
+   30 CONTINUE
+      KEV(I)=0
+   40 CONTINUE
+      DO 50 I=1,NFISS
+      KFISS(I)=KEV(KFISS(I))
+   50 CONTINUE
+*----
+*  FIRST LOOP OVER LOGY
+*----
+      DO 275 L=1,LOGY
+*----
+*  COMPUTE MATRICES A22**-1, A21, AND A12
+*----
+      DO 90 II=1,NSAT
+      I=KSAT(II)
+      IMAM1=0
+      IF(I.GT.1) IMAM1=IMA(I-1)
+      DO 60 JJ=1,NSAT
+      J=KSAT(JJ)
+      IF((J.LE.I).AND.(J.GT.I+IMAM1-MU1(I))) THEN
+         A22(II,JJ)=ADPL(MU1(I)-I+J,L)
+      ELSE IF((I.LE.J).AND.(I.GE.J-IMA(J)+MU1(J))) THEN
+         A22(II,JJ)=ADPL(MU1(J)+J-I,L)
+      ELSE
+         A22(II,JJ)=0.0
+      ENDIF
+   60 CONTINUE
+      JMAM1=0
+      DO 75 J=1,NVAR
+      J0=KEV(J)
+      IF(J0.EQ.0) GO TO 70
+      IF((J.LE.I).AND.(J.GT.I+IMAM1-MU1(I))) THEN
+         A21(II,J0,L)=ADPL(MU1(I)-I+J,L)
+      ELSE IF((I.LE.J).AND.(I.GE.J-IMA(J)+MU1(J))) THEN
+         A21(II,J0,L)=ADPL(MU1(J)+J-I,L)
+      ELSE
+         A21(II,J0,L)=0.0
+      ENDIF
+      IF((I.LE.J).AND.(I.GT.J+JMAM1-MU1(J))) THEN
+         A12(J0,II,L)=ADPL(MU1(J)-J+I,L)
+      ELSE IF((J.LE.I).AND.(J.GE.I-IMA(I)+MU1(I))) THEN
+         A12(J0,II,L)=ADPL(MU1(I)+I-J,L)
+      ELSE
+         A12(J0,II,L)=0.0
+      ENDIF
+   70 JMAM1=IMA(J)
+   75 CONTINUE
+      IF(I.GT.NSUPL) THEN
+         DO 80 K=1,NFISS
+         A21(II,KFISS(K),L)=A21(II,KFISS(K),L)+YSF(K,I-NSUPL,L)
+   80    CONTINUE
+      ENDIF
+   90 CONTINUE
+      CALL ALINV(NSAT,A22,NSAT,IER)
+      IF(IER.NE.0) CALL XABORT('EVOSAT: SINGULAR MATRIX.')
+*----
+*  COMPUTE VECTOR YSTG ANT YSAT
+*----
+      IF(L.EQ.1) THEN
+*        BEGINNING-OF-STAGE DIRAC DELTA CONTRIBUTIONS:
+         DO 100 I=1,NSAT
+         YSAT(I)=YST1(KSAT(I))
+  100    CONTINUE
+         DO 110 I=1,NVAR
+         IF(KEV(I).GT.0) YSTG(KEV(I))=YST1(I)
+  110    CONTINUE
+         IF(IDIRAC.EQ.0) THEN
+            DO 125 I=1,NSAT
+            GAR(I)=BDPL(KSAT(I),L)
+            DO 120 J=1,NVAR-NSAT
+            GAR(I)=GAR(I)+A21(I,J,L)*YSTG(J)
+  120       CONTINUE
+  125       CONTINUE
+            DO 135 I=1,NSAT
+            YSAT(I)=0.0
+            DO 130 J=1,NSAT
+            YSAT(I)=YSAT(I)-A22(I,J)*GAR(J)
+  130       CONTINUE
+  135       CONTINUE
+            GO TO 220
+         ENDIF
+         ITER=0
+  140    ITER=ITER+1
+         IF(ITER.GT.50) CALL XABORT('EVOSAT: CONVERGENCE FAILURE.')
+         DO 155 I=1,NSAT
+         GAR(I)=BDPL(KSAT(I),L)
+         DO 150 J=1,NVAR-NSAT
+         GAR(I)=GAR(I)+A21(I,J,L)*YSTG(J)
+  150    CONTINUE
+  155    CONTINUE
+         ERR1=0.0
+         ERR2=0.0
+         DO 170 I=1,NSAT
+         ZCOMP=YSAT(I)
+         YSAT(I)=0.0
+         DO 160 J=1,NSAT
+         YSAT(I)=YSAT(I)-A22(I,J)*GAR(J)
+  160    CONTINUE
+         ERR1=MAX(ERR1,ABS(ZCOMP-YSAT(I)))
+         ERR2=MAX(ERR2,ABS(YSAT(I)))
+  170    CONTINUE
+         DO 185 I=1,NSAT
+         GAR(I)=0.0
+         DO 180 J=1,NSAT
+         GAR(I)=GAR(I)-A22(I,J)*(YST1(KSAT(J))-YSAT(J))
+  180    CONTINUE
+  185    CONTINUE
+         DO 190 I=1,NVAR
+         IF(KEV(I).GT.0) YSTG(KEV(I))=YST1(I)
+  190    CONTINUE
+         DO 210 I=1,NVAR-NSAT
+         DO 200 J=1,NSAT
+         YSTG(I)=YSTG(I)+A12(I,J,L)*GAR(J)
+  200    CONTINUE
+         ERR2=MAX(ERR2,ABS(YSTG(I)))
+  210    CONTINUE
+         IF(ERR1.LE.EPS*ERR2) GO TO 220
+         GO TO 140
+      ENDIF
+*----
+*  COMPUTE MATRICES A21 AND BGAR
+*----
+  220 DO 235 I=1,NSAT
+      GAR(I)=0.0
+      DO 230 J=1,NSAT
+      GAR(I)=GAR(I)-A22(I,J)*BDPL(KSAT(J),L)
+  230 CONTINUE
+  235 CONTINUE
+      BGAR(:NVAR-NSAT)=0.0
+      DO 240 I=1,NVAR
+      IF(KEV(I).GT.0) BGAR(KEV(I))=BDPL(I,L)
+  240 CONTINUE
+      DO 255 I=1,NVAR-NSAT
+      DO 250 J=1,NSAT
+      BGAR(I)=BGAR(I)+A12(I,J,L)*GAR(J)
+  250 CONTINUE
+  255 CONTINUE
+      DO 272 J=1,NVAR-NSAT
+      BDPL(J,L)=BGAR(J)
+      IF(L.EQ.1) YST1(J)=YSTG(J)
+      DO 260 K=1,NSAT
+      GAR(K)=A21(K,J,L)
+  260 CONTINUE
+      DO 271 I=1,NSAT
+      A21(I,J,L)=0.0
+      DO 270 K=1,NSAT
+      A21(I,J,L)=A21(I,J,L)+A22(I,K)*GAR(K)
+  270 CONTINUE
+  271 CONTINUE
+  272 CONTINUE
+*
+  275 CONTINUE
+*----
+*  DETERMINE THE PROFILE PATTERN OF THE LUMPED DEPLETION MATRIX.
+*----
+      NSUPLG=NSUPL
+      DO 280 I=1,NVAR
+      IF((KEV(I).EQ.0).AND.(I.LE.NSUPL)) NSUPLG=NSUPLG-1
+  280 CONTINUE
+      NSUPFG=NVAR-NSAT-NSUPLG
+      MGAR(:NVAR-NSAT)=1
+      IGAR(:NVAR-NSAT)=1
+      IMAM1=0
+      DO 305 I=1,NVAR
+      IKEV=KEV(I)
+      IF(IKEV.EQ.0) GO TO 300
+      DO 290 J=1,NVAR
+      JKEV=KEV(J)
+      IF(JKEV.EQ.0) GO TO 290
+      IF((J.LE.I).AND.(J.GT.I+IMAM1-MU1(I))) THEN
+         MGAR(IKEV)=MAX(MGAR(IKEV),IKEV-JKEV+1)
+      ELSE IF((I.LE.J).AND.(I.GE.J-IMA(J)+MU1(J))) THEN
+         IGAR(JKEV)=MAX(IGAR(JKEV),JKEV-IKEV+1)
+      ENDIF
+  290 CONTINUE
+  300 IMAM1=IMA(I)
+  305 CONTINUE
+      DO 335 J=1,NVAR-NSAT
+      JIFI=0
+      DO 310 IFI=1,NFISS
+      IF(J.EQ.KFISS(IFI)) JIFI=IFI
+  310 CONTINUE
+      DO 330 I=1,NVAR-NSAT
+      IF((I.GT.NSUPLG).AND.(JIFI.GT.0)) GO TO 330
+      LTEST=.FALSE.
+      DO 325 L=1,LOGY
+      DO 320 K=1,NSAT
+      LTEST=LTEST.OR.(A12(I,K,L)*A21(K,J,L).NE.0.0)
+  320 CONTINUE
+  325 CONTINUE
+      IF(LTEST.AND.(J.LE.I)) THEN
+         MGAR(I)=MAX(MGAR(I),I-J+1)
+      ELSE IF(LTEST) THEN
+         IGAR(J)=MAX(IGAR(J),J-I+1)
+      ENDIF
+  330 CONTINUE
+  335 CONTINUE
+      II=0
+      DO 340 I=1,NVAR-NSAT
+      II=II+MGAR(I)
+      MGAR(I)=II
+      II=II+IGAR(I)-1
+      IGAR(I)=II
+  340 CONTINUE
+      IF(IMPX.GT.8) WRITE(6,'(/27H EVOSAT: REAL SIZE OF ADPL=,I9,3H AL,
+     1 13HLOCATED SIZE=,I9,1H.)') IGAR(NVAR-NSAT),MAXA
+      IF(IGAR(NVAR-NSAT).GT.MAXA) THEN
+         WRITE(HSMG,'(24HEVOSAT: IGAR(NVAR-NSAT)=,I6,6H MAXA=,I6)')
+     1   IGAR(NVAR-NSAT),MAXA
+         CALL XABORT(HSMG)
+      ENDIF
+*----
+*  SECOND LOOP OVER LOGY
+*----
+      DO 540 L=1,LOGY
+*----
+*  COMPUTE MATRIX AGAR AND YIELDS YSFG.
+*----
+      AGAR(:IGAR(NVAR-NSAT))=0.0
+      IMAM1=0
+      DO 445 I=1,NVAR
+      IKEV=KEV(I)
+      IF(IKEV.EQ.0) GO TO 440
+      DO 420 J=1,NVAR
+      JKEV=KEV(J)
+      IF(JKEV.EQ.0) GO TO 420
+      IF((J.LE.I).AND.(J.GT.I+IMAM1-MU1(I))) THEN
+         AGAR(MGAR(IKEV)-IKEV+JKEV)=ADPL(MU1(I)-I+J,L)
+      ELSE IF((I.LE.J).AND.(I.GE.J-IMA(J)+MU1(J))) THEN
+         AGAR(MGAR(JKEV)+JKEV-IKEV)=ADPL(MU1(J)+J-I,L)
+      ENDIF
+  420 CONTINUE
+      IF(I.GT.NSUPL) THEN
+         DO 430 K=1,NFISS
+         YSFG(K,IKEV-NSUPLG)=YSF(K,I-NSUPL,L)
+  430    CONTINUE
+      ENDIF
+  440 IMAM1=IMA(I)
+  445 CONTINUE
+      DO 495 J=1,NVAR-NSAT
+      JIFI=0
+      DO 450 IFI=1,NFISS
+      IF(J.EQ.KFISS(IFI)) JIFI=IFI
+  450 CONTINUE
+      IMAM1=0
+      DO 490 I=1,NVAR-NSAT
+      IF((I.GT.NSUPLG).AND.(JIFI.GT.0)) GO TO 480
+      IF((J.LE.I).AND.(J.GT.I+IMAM1-MGAR(I))) THEN
+         DO 460 K=1,NSAT
+         AGAR(MGAR(I)-I+J)=AGAR(MGAR(I)-I+J)-A12(I,K,L)*A21(K,J,L)
+  460    CONTINUE
+      ELSE IF((I.LE.J).AND.(I.GE.J-IGAR(J)+MGAR(J))) THEN
+         DO 470 K=1,NSAT
+         AGAR(MGAR(J)+J-I)=AGAR(MGAR(J)+J-I)-A12(I,K,L)*A21(K,J,L)
+  470    CONTINUE
+      ENDIF
+  480 IMAM1=IGAR(I)
+  490 CONTINUE
+  495 CONTINUE
+      DO 510 I=NSUPLG+1,NVAR-NSAT
+      DO 505 IFI=1,NFISS
+      J=KFISS(IFI)
+      DO 500 K=1,NSAT
+      YSFG(IFI,I-NSUPLG)=YSFG(IFI,I-NSUPLG)-A12(I,K,L)*A21(K,J,L)
+  500 CONTINUE
+  505 CONTINUE
+  510 CONTINUE
+*----
+*  REPLACE THE ORIGINAL INFORMATION WITH THE LUMPED ONE
+*----
+      DO 520 I=1,IGAR(NVAR-NSAT)
+      ADPL(I,L)=AGAR(I)
+  520 CONTINUE
+      DO 535 I=1,NFISS
+      DO 530 J=1,NSUPFG
+      YSF(I,J,L)=YSFG(I,J)
+  530 CONTINUE
+  535 CONTINUE
+  540 CONTINUE
+      DO 550 I=1,NVAR-NSAT
+      IMA(I)=IGAR(I)
+      MU1(I)=MGAR(I)
+  550 CONTINUE
+      RETURN
+      END