Initial commit from Polytechnique Montreal

1 files changed, 229 insertions, 0 deletions

diff --git a/Dragon/src/EVOSOL.f b/Dragon/src/EVOSOL.f
new file mode 100644
index 0000000..d256faa
--- /dev/null
+++ b/Dragon/src/EVOSOL.f
@@ -0,0 +1,229 @@
+*DECK EVOSOL
+      SUBROUTINE EVOSOL(IMPX,LCOOL,NVAR,NREAC,NDFP,NPAR,NFISS,XT,EPS1,
+     1 EXPMAX,H1,ITYPE,IDIRAC,DCR,KPAR,BPAR,KFISS,KPF,YIELD,LP,IEVOLB,
+     2 SIG1,SIG2,NVAR2,NFISS2,NSUPF2,MU1,IMA,MAXA,YDPL,ICHAIN)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Put the depletion matrix system in sparse storage mode for a single
+* depleting mixture. Solve this system between times XT(1) and XT(2).
+*
+*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).
+* LCOOL   out-of-core depletion flag (LCOOL=.true. to set flag).
+* NVAR    number of depleting nuclides.
+* NREAC   one plus the number of neutron-induced depletion reactions.
+* NDFP    number of direct fission products (fission fragments).
+* NPAR    maximum number of parent nuclides in the depletion chain.
+* NFISS   number of fissile isotopes producing fission products.
+* XT      initial and final time (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.
+* 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).
+* DCR     sum of radioactive decay constants in 10**-8/s
+* KPAR    position in chain of the parent nuclide and type of
+*         reaction.
+* BPAR    branching ratio for neutron induced reactions.
+* KFISS   position in chain of the fissile isotopes.
+* KPF     position in chain of the direct fission products (fission
+*         fragments).
+* YIELD   fission yields.
+* LP      index vector used to remove unused isotopes from the
+*         depletion chain.
+* IEVOLB  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.
+* SIG1    initial reaction rates for nuclide I:
+*         SIG1(I,1) fission reaction rate;
+*         SIG1(I,2) gamma reaction rate;
+*         SIG1(I,3) N2N reaction rate;
+*         ...;
+*         SIG1(I,NREAC) neutron-induced energy;
+*         SIG1(I,NREAC+1) decay energy released.
+* SIG2    final reaction rates.
+* NVAR2   number of used isotopes in the depletion chain, where
+*         NVAR2=max(LP(I)).
+* NFISS2  number of used fissile isotopes producing fission products.
+* NSUPF2  number of used fission products.
+* 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.
+* YDPL    initial/final number density of each depleting isotope.
+*         YDPL(NVAR+1,2) is the stage burnup increment.
+* ICHAIN  name of the used isotopes in the depletion chain.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      LOGICAL LCOOL
+      INTEGER IMPX,NVAR,NREAC,NDFP,NPAR,NFISS,ITYPE,IDIRAC,
+     1 KPAR(NPAR,NVAR),KFISS(NFISS),KPF(NDFP),LP(NVAR),IEVOLB(NVAR),
+     2 NVAR2,NFISS2,NSUPF2,MU1(NVAR2+1),IMA(NVAR2+1),MAXA,
+     3 ICHAIN(2,NVAR2+1)
+      REAL XT(2),EPS1,EXPMAX,H1,DCR(NVAR),BPAR(NPAR,NVAR),
+     1 YIELD(NFISS,NDFP),SIG1(NVAR+1,NREAC+1),SIG2(NVAR+1,NREAC+1),
+     2 YDPL(NVAR+1,2)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: LQ,KFISS2,IEVOL2
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: ADPL,BDPL,YDPL2
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: YSF
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(LQ(NFISS),KFISS2(NFISS2),IEVOL2(NVAR2+1))
+      ALLOCATE(ADPL(MAXA,2),BDPL(NVAR2+1,2),YSF(NFISS2,NSUPF2+1,2),
+     1 YDPL2(NVAR2+1,2))
+*----
+*  COMPUTE LQ AND KFISS2
+*----
+      I0=0
+      DO 10 I=1,NFISS
+      LQ(I)=0
+      IF(KFISS(I).EQ.0) GO TO 10
+      IF((LP(KFISS(I)).EQ.0).OR.(LP(KFISS(I)).GT.NVAR2)) GO TO 10
+      I0=I0+1
+      IF(I0.GT.NFISS2) CALL XABORT('EVOSOL: NFISS2 TOO SMALL.')
+      KFISS2(I0)=LP(KFISS(I))
+      LQ(I)=I0
+   10 CONTINUE
+*----
+*  BUILD THE SPARSE DEPLETION MATRICES ADPL, BDPL AND YSF
+*----
+      DO 80 IP=1,2
+      DO 20 I=1,IMA(NVAR2+1)
+      ADPL(I,IP)=0.0
+   20 CONTINUE
+      BDPL(NVAR2+1,IP)=0.0
+      DO 40 IS=1,NVAR
+      IF((LP(IS).EQ.0).OR.(LP(IS).GT.NVAR2)) GO TO 40
+      BDPL(LP(IS),IP)=0.0
+      SIGE=0.0
+      IF((.NOT.LCOOL).AND.(IP.EQ.1)) THEN
+         DO 25 IREAC=1,NREAC-1
+         SIGE=SIGE+SIG1(IS,IREAC)
+   25    CONTINUE
+      ELSE IF((.NOT.LCOOL).AND.(IP.EQ.2)) THEN
+         DO 26 IREAC=1,NREAC-1
+         SIGE=SIGE+SIG2(IS,IREAC)
+   26    CONTINUE
+      ENDIF
+      ADPL(MU1(LP(IS)),IP)=-SIGE-DCR(IS)
+      DO 30 IPAR=1,NPAR
+      IF(KPAR(IPAR,IS).EQ.0) GO TO 40
+      IF(KPAR(IPAR,IS).EQ.2) GO TO 30
+      JS=KPAR(IPAR,IS)/100
+      KT=KPAR(IPAR,IS)-JS*100
+      IF((LCOOL.AND.(KT.GE.2)).OR.(LP(JS).EQ.0).OR.(LP(JS).GT.NVAR2))
+     1 GO TO 30
+      SIGE=0.0
+      IF(KT.EQ.1) THEN
+         SIGE=BPAR(IPAR,IS)*DCR(JS)
+      ELSE IF((KT.GE.3).AND.(IP.EQ.1)) THEN
+         SIGE=BPAR(IPAR,IS)*SIG1(JS,KT-1)
+      ELSE IF((KT.GE.3).AND.(IP.EQ.2)) THEN
+         SIGE=BPAR(IPAR,IS)*SIG2(JS,KT-1)
+      ELSE
+         CALL XABORT('EVOSOL: UNKNOWN REACTION.')
+      ENDIF
+      IF(LP(JS).LE.LP(IS)) THEN
+         ADPL(MU1(LP(IS))-LP(IS)+LP(JS),IP)=
+     1   ADPL(MU1(LP(IS))-LP(IS)+LP(JS),IP)+SIGE
+      ELSE
+         ADPL(MU1(LP(JS))+LP(JS)-LP(IS),IP)=
+     1   ADPL(MU1(LP(JS))+LP(JS)-LP(IS),IP)+SIGE
+      ENDIF
+   30 CONTINUE
+   40 CONTINUE
+      IF(LCOOL) THEN
+         DO 51 IFIS=1,NFISS2
+         DO 50 ISUPF=1,NSUPF2+1
+         YSF(IFIS,ISUPF,IP)=0.0
+   50    CONTINUE
+   51    CONTINUE
+      ELSE
+*        ADD ONE EQUATION TO COMPUTE THE BURNUP.
+         DO 55 JS=1,NVAR
+         IF((LP(JS).EQ.0).OR.(LP(JS).GT.NVAR2)) GO TO 55
+         IF(IP.EQ.1) THEN
+            ADPL(MU1(NVAR2+1)-(NVAR2+1)+LP(JS),IP)=SIG1(JS,NREAC)+
+     &      SIG1(JS,NREAC+1)
+         ELSE IF(IP.EQ.2) THEN
+            ADPL(MU1(NVAR2+1)-(NVAR2+1)+LP(JS),IP)=SIG2(JS,NREAC)+
+     &      SIG2(JS,NREAC+1)
+         ENDIF
+   55    CONTINUE
+*
+*        ADD THE FISSION YIELD CONTRIBUTIONS.
+         DO 70 IFIS=1,NFISS
+         IF(LQ(IFIS).EQ.0) GO TO 70
+         DO 56 ISUPF=1,NSUPF2+1
+         YSF(LQ(IFIS),ISUPF,IP)=0.0
+   56    CONTINUE
+         IF(NSUPF2.EQ.0) GO TO 70
+         DO 60 ISUPF=1,NDFP
+         LPP=LP(KPF(ISUPF))
+         IF((LPP.EQ.0).OR.(LPP.GT.NVAR2)) GO TO 60
+         IF(LPP+NSUPF2.LE.NVAR2) CALL XABORT('EVOSOL: FAILURE.')
+         IF(IP.EQ.1) THEN
+            YSF(LQ(IFIS),LPP-NVAR2+NSUPF2,IP)=YIELD(IFIS,ISUPF)*
+     1      SIG1(KFISS(IFIS),1)
+         ELSE IF(IP.EQ.2) THEN
+            YSF(LQ(IFIS),LPP-NVAR2+NSUPF2,IP)=YIELD(IFIS,ISUPF)*
+     1      SIG2(KFISS(IFIS),1)
+         ENDIF
+   60    CONTINUE
+   70    CONTINUE
+      ENDIF
+   80 CONTINUE
+*----
+*  SOLVE THE DEPLETION SYSTEM. EQUATION NVAR2+1 IS USED TO COMPUTE
+*  THE BURNUP
+*----
+      YDPL2(:NVAR2+1,1)=0.0
+      IEVOL2(:NVAR2+1)=0
+      DO 90 IS=1,NVAR
+      IF((LP(IS).EQ.0).OR.(LP(IS).GT.NVAR2)) GO TO 90
+      YDPL2(LP(IS),1)=YDPL(IS,1)
+      IEVOL2(LP(IS))=IEVOLB(IS)
+   90 CONTINUE
+      CALL EVODPL(IMPX,YDPL2(1,1),NVAR2+1,XT,EPS1,EXPMAX,H1,ITYPE,
+     1 IDIRAC,IEVOL2,MU1,IMA,MAXA,NSUPF2+1,NFISS2,KFISS2,YSF,ADPL,
+     2 BDPL,ICHAIN)
+      YDPL(NVAR+1,2)=YDPL2(NVAR2+1,2)
+      DO 100 IS=NVAR,1,-1
+      IF((LP(IS).EQ.0).OR.(LP(IS).GT.NVAR2)) THEN
+         YDPL(IS,2)=YDPL(IS,1)
+      ELSE
+         YDPL(IS,1)=YDPL2(LP(IS),1)
+         YDPL(IS,2)=MAX(0.0,YDPL2(LP(IS),2))
+      ENDIF
+  100 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(YDPL2,YSF,BDPL,ADPL)
+      DEALLOCATE(IEVOL2,KFISS2,LQ)
+      RETURN
+      END