Initial commit from Polytechnique Montreal

1 files changed, 425 insertions, 0 deletions

diff --git a/Dragon/src/EVODPL.f b/Dragon/src/EVODPL.f
new file mode 100644
index 0000000..7868b98
--- /dev/null
+++ b/Dragon/src/EVODPL.f
@@ -0,0 +1,425 @@
+*DECK EVODPL
+      SUBROUTINE EVODPL(IMPX,YDPL,NVAR,XT,EPS1,EXPMAX,H1,ITYPE,IDIRAC,
+     1 IEVOL2,MU1,IMA,MAXA,NSUPF,NFISS,KFISS,YSF,ADPL,BDPL,ICHAIN)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Multi-purpose driver for solving the isotopic depletion equations,
+* taking into account the saturation phenomena.
+*
+*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/output
+* IMPX    print flag (equal to zero for no print).
+* YDPL    initial/final number densities.
+* NVAR    number of nuclides in the complete depletion chain.
+* XT      initial and final value of the independent variable.
+* EPS1    required accuracy for the ODE solver.
+* EXPMAX  saturation limit. A nuclide is saturating if
+*         -ADPL(MU1(I))*(XT(2)-XT(1)).GT.EXPMAX. Suggested value:
+*         EXPMAX=80.0. EXPMAX=0.0 means that the saturation model is
+*         not used.
+* H1      guessed first stepsize.
+* ITYPE   type of ODE solution:
+*         =1 fifth-order Runge-Kutta method;
+*         =2 fourth-order Kaps-Rentrop method.
+* IDIRAC  saturation model flag (=1 to use Dirac function contributions
+*         in the saturating nuclide number densities.
+* IEVOL2  flag making an isotope non-depleting:
+*         =1 to force an isotope to be non-depleting;
+*         =2 to force an isotope to be depleting;
+*         =3 to force an isotope to be at saturation.
+* MU1     position of each diagonal element in matrix ADPL.
+* IMA     position of the first non-zero column element in matrix ADPL.
+* MAXA    first dimension of matrix ADPL.
+* NSUPF   number of depleting fission products.
+* NFISS   number of fissile isotopes producing fission products.
+* KFISS   position in chain of the fissile isotopes.
+* YSF     initial/final product of the fission yields and fission
+*         rates.
+* ADPL    initial/final depletion matrix.
+* BDPL    initial/final depletion source.
+* ICHAIN  name of the isotopes in the depletion chain.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IMPX,NVAR,ITYPE,IDIRAC,IEVOL2(NVAR),MU1(NVAR),IMA(NVAR),
+     1 MAXA,NSUPF,NFISS,KFISS(NFISS),ICHAIN(2,NVAR)
+      REAL YDPL(NVAR,2),XT(2),EPS1,EXPMAX,H1,YSF(NFISS,NSUPF,2),
+     1 ADPL(MAXA,2),BDPL(NVAR,2)
+*----
+*  LOCAL VARIABLES
+*----
+      LOGICAL LSAT
+      CHARACTER*2 SHOW(120,120)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: KSAT,IPERM,MU12,IMA2,KFIS2
+      REAL, ALLOCATABLE, DIMENSION(:) :: YST1,YSAT
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: ADPL2,BDPL2,BDPL3
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: YSF2,YSF3
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(KSAT(NVAR),IPERM(NVAR),MU12(NVAR),IMA2(NVAR),
+     1 KFIS2(NFISS))
+      ALLOCATE(YST1(NVAR),YSAT(NVAR),ADPL2(MAXA,2),BDPL2(NVAR,2),
+     1 YSF2(NFISS,NSUPF,2))
+*----
+*  COMPUTE THE LUMPING INDEX VECTOR IPERM
+*----
+      DO 10 I=1,NVAR
+      IPERM(I)=I
+   10 CONTINUE
+      NTER=0
+   20 NTER=NTER+1
+      INDSAT=0
+      DO 30 I=1,NVAR
+      IF(IPERM(I).GE.0) THEN
+         LSAT=(IEVOL2(I).EQ.3).AND.(EXPMAX.GT.0.0)
+         IF(EXPMAX.GT.0.0) THEN
+            LSAT=LSAT.OR.((ABS(ADPL(MU1(I),1)*(XT(2)-XT(1))).GT.EXPMAX)
+     >      .AND.(ABS(ADPL(MU1(I),2)*(XT(2)-XT(1))).GT.EXPMAX))
+         ENDIF
+         IF(LSAT) THEN
+            DO 25 II=1,NFISS
+            IF(I.EQ.KFISS(II)) GO TO 30
+   25       CONTINUE
+            IPERM(I)=0
+            IF(INDSAT.EQ.0) THEN
+               IF(IMPX.GT.5) WRITE(6,'(17H EVODPL: ISOTOPE ,2A4,
+     1         18H IS SATURATING(1).)') ICHAIN(1,I),ICHAIN(2,I)
+               IPERM(I)=-NTER
+               INDSAT=I
+            ENDIF
+         ENDIF
+      ENDIF
+   30 CONTINUE
+      IF(INDSAT.EQ.0) GO TO 60
+      DO 50 I=INDSAT+1,NVAR
+      JMN=I-MU1(I)+IMA(I-1)+1
+      IMN=I-IMA(I)+MU1(I)
+      DO 40 J=MIN(JMN,IMN),I-1
+      IF((IPERM(I).EQ.0).AND.(IPERM(J).EQ.-NTER)) THEN
+         IF(IMPX.GT.5) WRITE(6,'(17H EVODPL: ISOTOPE ,2A4,
+     1   18H IS SATURATING(2).)') ICHAIN(1,I),ICHAIN(2,I)
+         IPERM(I)=-NTER
+         GO TO 50
+      ENDIF
+   40 CONTINUE
+   50 CONTINUE
+      GO TO 20
+   60 NTER=NTER-1
+      N=0
+      DO 70 I=1,NVAR
+      IF(IPERM(I).GT.0) THEN
+         N=N+1
+         IPERM(I)=N
+      ENDIF
+   70 CONTINUE
+      IF(IMPX.GT.3) THEN
+         WRITE(6,400) NVAR,XT(1),XT(2),EPS1,H1,ITYPE,NTER,NVAR-N,NFISS,
+     1   NSUPF,(IPERM(I),I=1,NVAR)
+         WRITE(6,410) (YDPL(I,1),I=1,NVAR)
+      ENDIF
+      IF(IMPX.GT.5) THEN
+         NVARM=MIN(NVAR,120)
+         WRITE (6,'(//34H EVODPL: DEPLETION MATRIX PROFILE:/)')
+         DO 85 I=1,NVARM
+         DO 80 J=1,NVARM
+         SHOW(I,J)=' '
+   80    CONTINUE
+   85    CONTINUE
+         IMAM1=0
+         DO 120 I=1,NVARM
+         DO 90 J=I-MU1(I)+IMAM1+1,I-1
+         SHOW(I,J)='*'
+   90    CONTINUE
+         DO 100 J=I-IMA(I)+MU1(I),I-1
+         SHOW(J,I)='*'
+  100    CONTINUE
+         IF(I.GT.NVAR-NSUPF) THEN
+            DO 110 K=1,NFISS
+            KFI=KFISS(K)
+            IF((KFI.GT.0).AND.(KFI.LE.120)) SHOW(I,KFI)='-'
+  110       CONTINUE
+         ENDIF
+         SHOW(I,I)='+'
+         IMAM1=IMA(I)
+  120    CONTINUE
+         DO 130 I=1,NVARM
+         WRITE (6,'(1X,I4,1X,2A4,1X,120A2)') I,ICHAIN(1,I),ICHAIN(2,I),
+     1   (SHOW(I,J),J=1,NVARM)
+  130    CONTINUE
+         IF(NVAR.GT.120) 
+     >   WRITE(6,'(34H MATRIX TRUNCATED TO 120 ELEMENTS.)')
+         IF(IMPX .GE. 1000) THEN
+           write(6,'(A)') 'ORIGINAL DEPLETION SYSTEM'
+           write(6,'(3I10)') NVAR,NFISS,NSUPF
+           write(6,'(A6)') 'MU1   '
+           write(6,'(20I5)') (MU1(I),I=1,NVAR)
+           write(6,'(A6)') 'IMA   '
+           write(6,'(20I5)') (IMA(I),I=1,NVAR)
+           write(6,'(A6)') 'ADPL1 '
+           write(6,'(1P,5E20.12)') (ADPL(I,1),I=1,IMA(NVAR))
+           write(6,'(A6)') 'BDPL1 '
+           write(6,'(1P,5E20.12)') (BDPL(I,1),I=1,NVAR)
+           write(6,'(A6)') 'KFISS '
+           write(6,'(20I5)') (KFISS(K),K=1,NFISS)
+           write(6,'(A6)') 'YSF1  '
+           write(6,'(1P,5E20.12)') ((YSF(I,J,1),I=1,NFISS),J=1,NSUPF)
+         ENDIF
+      ENDIF
+*----
+*  LUMPING OF THE DEPLETION MATRICES
+*----
+      DO 135 IFI=1,NFISS
+      KFIS2(IFI)=KFISS(IFI)
+  135 CONTINUE
+      DO 140 I=1,NVAR
+      YST1(I)=YDPL(I,1)
+      MU12(I)=MU1(I)
+      IMA2(I)=IMA(I)
+  140 CONTINUE
+      DO 162 L=1,2
+      DO 145 I=1,NVAR
+      BDPL2(I,L)=BDPL(I,L)
+  145 CONTINUE
+      DO 150 I=1,IMA(NVAR)
+      ADPL2(I,L)=ADPL(I,L)
+  150 CONTINUE
+      DO 161 I=1,NFISS
+      DO 160 J=1,NSUPF
+      YSF2(I,J,L)=YSF(I,J,L)
+  160 CONTINUE
+  161 CONTINUE
+  162 CONTINUE
+      NVAR2=NVAR
+      NSUPF2=NSUPF
+      DO 180 ITER=1,NTER
+      I0=0
+      NSAT=0
+      DO 170 I=1,NVAR
+      IF((IPERM(I).GT.0).OR.(IPERM(I).LT.-ITER)) THEN
+        I0=I0+1
+      ELSE IF(IPERM(I).EQ.-ITER) THEN
+        I0=I0+1
+        NSAT=NSAT+1
+        KSAT(NSAT)=I0
+      ENDIF
+  170 CONTINUE
+      IF(I0.NE.NVAR2) CALL XABORT('EVODPL: ALGORITHM FAILURE 1.')
+      MAXB=NVAR
+      MAXY=NSUPF
+      CALL EVOSAT(IMPX,MAXA,MAXB,MAXY,2,NSAT,NVAR2,KSAT,YST1,YSAT,MU12,
+     1 IMA2,NSUPF2,NFISS,IDIRAC,KFIS2,YSF2(1,1,1),ADPL2(1,1),BDPL2(1,1),
+     2 NSUPF3)
+      NVAR2=NVAR2-NSAT
+      NSUPF2=NSUPF3
+      NSAT=0
+      I0=0
+      DO 175 I=1,NVAR
+      IF((IPERM(I).GT.0).OR.(IPERM(I).LT.-ITER)) THEN
+         I0=I0+1
+         YDPL(I,1)=YST1(I0)
+      ELSE IF(IPERM(I).EQ.-ITER) THEN
+         NSAT=NSAT+1
+         YDPL(I,1)=YSAT(NSAT)
+      ENDIF
+  175 CONTINUE
+  180 CONTINUE
+      IF(IMPX.GT.4) WRITE(6,420) (YDPL(I,1),I=1,NVAR)
+*----
+*  SOLUTION OF THE LUMPED DEPLETION SYSTEM
+*----
+      DO 185 I=1,NVAR
+      YDPL(I,2)=YDPL(I,1)
+  185 CONTINUE
+      IF(NVAR2.EQ.0) GO TO 315
+      DO 190 I=1,NVAR2
+      FACT=(BDPL2(I,2)-BDPL2(I,1))/(XT(2)-XT(1))
+      BDPL2(I,1)=BDPL2(I,1)-FACT*XT(1)
+      BDPL2(I,2)=FACT
+  190 CONTINUE
+      DO 200 I=1,IMA2(NVAR2)
+      FACT=(ADPL2(I,2)-ADPL2(I,1))/(XT(2)-XT(1))
+      ADPL2(I,1)=ADPL2(I,1)-FACT*XT(1)
+      ADPL2(I,2)=FACT
+  200 CONTINUE
+      DO 215 I=1,NFISS
+      DO 210 J=1,NSUPF2
+      FACT=(YSF2(I,J,2)-YSF2(I,J,1))/(XT(2)-XT(1))
+      YSF2(I,J,1)=YSF2(I,J,1)-FACT*XT(1)
+      YSF2(I,J,2)=FACT
+  210 CONTINUE
+  215 CONTINUE
+      IF(IMPX.GT.4) THEN
+         WRITE(6,430) NSUPF2
+         WRITE(6,440) (YST1(I),I=1,NVAR2)
+      ENDIF
+      IF(IMPX.GT.5) THEN
+         NVARM=MIN(NVAR2,120)
+         WRITE (6,'(//41H EVODPL: LUMPED DEPLETION MATRIX PROFILE:/)')
+         DO 225 I=1,NVARM
+         DO 220 J=1,NVARM
+         SHOW(I,J)=' '
+  220    CONTINUE
+  225    CONTINUE
+         IMAM1=0
+         DO 260 I=1,NVARM
+         DO 230 J=I-MU12(I)+IMAM1+1,I-1
+         SHOW(I,J)='*'
+  230    CONTINUE
+         DO 240 J=I-IMA2(I)+MU12(I),I-1
+         SHOW(J,I)='*'
+  240    CONTINUE
+         IF(I.GT.NVAR2-NSUPF2) THEN
+            DO 250 K=1,NFISS
+            KFI=KFIS2(K)
+            IF((KFI.GT.0).AND.(KFI.LE.60)) SHOW(I,KFI)='-'
+  250       CONTINUE
+         ENDIF
+         SHOW(I,I)='+'
+         IMAM1=IMA2(I)
+  260    CONTINUE
+         DO 270 I=1,NVARM
+         WRITE (6,'(1X,I4,1X,2A4,1X,120A2)') I,ICHAIN(1,I),ICHAIN(2,I),
+     1   (SHOW(I,J),J=1,NVARM)
+  270    CONTINUE
+         IF(NVAR.GT.120) 
+     >   WRITE(6,'(34H MATRIX TRUNCATED TO 120 ELEMENTS.)')
+         IF(IMPX .GE. 1000) THEN
+           write(6,'(A)') 'LUMPED DEPLETION SYSTEM'
+           write(6,'(3I10)') NVAR2,NFISS,NSUPF2
+           write(6,'(A6)') 'MU1   '
+           write(6,'(20I5)') (MU12(I),I=1,NVAR2)
+           write(6,'(A6)') 'IMA   '
+           write(6,'(20I5)') (IMA2(I),I=1,NVAR2)
+           write(6,'(A6)') 'ADPL2 '
+           write(6,'(1P,5E20.12)') (ADPL2(I,1),I=1,IMA2(NVAR2))
+           write(6,'(A6)') 'BDPL2 '
+           write(6,'(1P,5E20.12)') (BDPL2(I,1),I=1,NVAR2)
+           write(6,'(A6)') 'KFISS '
+           write(6,'(20I5)') (KFIS2(K),K=1,NFISS)
+           write(6,'(A6)') 'YSF1  '
+           write(6,'(1P,5E20.12)') ((YSF2(I,J,1),I=1,NFISS),J=1,NSUPF2)
+         ENDIF
+      ENDIF
+      ALLOCATE(BDPL3(NVAR2,2),YSF3(NFISS,NSUPF2,2))
+      DO 280 I=1,NVAR2
+      BDPL3(I,1)=BDPL2(I,1)
+      BDPL3(I,2)=BDPL2(I,2)
+  280 CONTINUE
+      DO 295 I=1,NFISS
+      DO 290 J=1,NSUPF2
+      YSF3(I,J,1)=YSF2(I,J,1)
+      YSF3(I,J,2)=YSF2(I,J,2)
+  290 CONTINUE
+  295 CONTINUE
+      CALL EVOODE(YST1,NVAR2,XT(1),XT(2),EPS1,H1,NOK,NBAD,ITYPE,MU12,
+     1 IMA2,MAXA,NSUPF2,NFISS,KFIS2,YSF3,ADPL2,BDPL3)
+      DEALLOCATE(YSF3,BDPL3)
+      IF(IMPX.GT.4) THEN
+         WRITE(6,450) (YST1(I),I=1,NVAR2)
+         IF(ITYPE.LE.2) WRITE(6,'(13H EVODPL: NOK=,I5,6H NBAD=,I5)')
+     1   NOK,NBAD
+      ENDIF
+      DO 310 I=1,NVAR
+      IF(IPERM(I).GT.0) YDPL(I,2)=YST1(IPERM(I))
+  310 CONTINUE
+*----
+*  COMPUTE NUMBER DENSITIES OF THE SATURATED ISOTOPES
+*----
+  315 IF(NTER.EQ.0) GO TO 370
+      DO 320 I=1,NVAR
+      YST1(I)=YDPL(I,2)
+      BDPL2(I,2)=BDPL(I,2)
+      MU12(I)=MU1(I)
+      IMA2(I)=IMA(I)
+  320 CONTINUE
+      DO 330 I=1,IMA(NVAR)
+      ADPL2(I,2)=ADPL(I,2)
+  330 CONTINUE
+      DO 345 I=1,NFISS
+      KFIS2(I)=KFISS(I)
+      DO 340 J=1,NSUPF
+      YSF2(I,J,2)=YSF(I,J,2)
+  340 CONTINUE
+  345 CONTINUE
+      NVAR2=NVAR
+      NSUPF2=NSUPF
+      DO 365 ITER=1,NTER
+      I0=0
+      NSAT=0
+      DO 350 I=1,NVAR
+      IF((IPERM(I).GT.0).OR.(IPERM(I).LT.-ITER)) THEN
+        I0=I0+1
+      ELSE IF(IPERM(I).EQ.-ITER) THEN
+        I0=I0+1
+        NSAT=NSAT+1
+        KSAT(NSAT)=I0
+      ENDIF
+  350 CONTINUE
+      IF(I0.NE.NVAR2) CALL XABORT('EVODPL: ALGORITHM FAILURE 2.')
+      MAXB=NVAR
+      MAXY=NSUPF
+      CALL EVOSAT(IMPX,MAXA,MAXB,MAXY,1,NSAT,NVAR2,KSAT,YST1,YSAT,MU12,
+     1 IMA2,NSUPF2,NFISS,IDIRAC,KFIS2,YSF2(1,1,2),ADPL2(1,2),BDPL2(1,2),
+     2 NSUPF3)
+      IF(IMPX.GT.4) WRITE(6,425) ITER,(YSAT(I),I=1,NSAT)
+      NVAR2=NVAR2-NSAT
+      NSUPF2=NSUPF3
+      NSAT=0
+      I0=0
+      DO 360 I=1,NVAR
+      IF((IPERM(I).GT.0).OR.(IPERM(I).LT.-ITER)) THEN
+        I0=I0+1
+        YDPL(I,2)=YST1(I0)
+      ELSE IF(IPERM(I).EQ.-ITER) THEN
+        NSAT=NSAT+1
+        YDPL(I,2)=YSAT(NSAT)
+      ENDIF
+  360 CONTINUE
+  365 CONTINUE
+  370 IF(IMPX.GT.3) WRITE(6,460) (YDPL(I,2),I=1,NVAR)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(YSF2,BDPL2,ADPL2,YSAT,YST1)
+      DEALLOCATE(KFIS2,IMA2,MU12,IPERM,KSAT)
+      RETURN
+*
+  400 FORMAT(//45H EVODPL: SOLUTION OF THE DEPLETION EQUATIONS.//14X,
+     1 25HTOTAL NUMBER OF NUCLIDES=,I5/26X,13HINITIAL TIME=,1P,E12.4/
+     2 28X,11HFINAL TIME=,E12.4/15X,24HACCURACY FOR ODE SOLVER=,E12.4/
+     3 16X,23HGUESSED FIRST STEPSIZE=,E12.4,0P/22X,17HTYPE OF SOLUTION=,
+     4 I3/39H NUMBER OF GROUP OF SATURATED NUCLIDES=,I5/10X,
+     5 29HNUMBER OF SATURATED NUCLIDES=,I5/12X,19HNUMBER OF FISSILE N,
+     6 8HUCLIDES=,I5/12X,27HNUMBER OF FISSION PRODUCTS=,I5//
+     7 22H LUMPING INDEX VECTOR:/(1X,20I5))
+  410 FORMAT(/48H EVODPL: INITIAL VALUES OF THE DEPLETION SYSTEM:/
+     1 (1X,1P,10E12.4))
+  420 FORMAT(/53H EVODPL: SATURATED INITIAL CONDITIONS OF THE DEPLETIO,
+     1 9HN SYSTEM:/(1X,1P,10E12.4))
+  425 FORMAT(/51H EVODPL: FINAL VALUES OF THE SATURATED NUCLIDES IN ,
+     1 9HGROUP NO.,I5//(1X,1P,10E12.4))
+  430 FORMAT(/42H NUMBER OF NON-SATURATED FISSION PRODUCTS=,I5)
+  440 FORMAT(/55H EVODPL: INITIAL VALUES OF THE LUMPED DEPLETION SYSTEM:
+     1 /(1X,1P,10E12.4))
+  450 FORMAT(/53H EVODPL: ODE SOLUTION OF THE LUMPED DEPLETION SYSTEM:/
+     1 (1X,1P,10E12.4))
+  460 FORMAT(/42H EVODPL: SOLUTION OF THE DEPLETION SYSTEM:/
+     1 (1X,1P,10E12.4))
+      END