diff options
Diffstat (limited to 'Dragon/src/EVOSOL.f')
| -rw-r--r-- | Dragon/src/EVOSOL.f | 229 |
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 |