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
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
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
|
