1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
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
|
