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
|
*DECK LIBWET
SUBROUTINE LIBWET(MAXR,NEL,NBESP,NSTATE,NMDEPL,ITNAM,ISTATE,
> MATNO,KPAX,BPAX)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Identify the depleting isotopes by type and reorder them (recompute
* the KPAX and BPAX matrices).
*
*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): G. Marleau
*
*Parameters: input
* MAXR number of reaction types.
* NEL number of isotopes on library.
* NBESP number of energy-dependent fission yield matrices.
* NSTATE number of parameters in the state vector.
* NMDEPL names of reactions:
* NMDEPL(1)='DECAY'; NMDEPL(2)='NFTOT';
* NMDEPL(3)='NG' ; NMDEPL(4)='N2N';
* etc.
* ITNAM reactive isotope names in chain.
*
*Parameters: output
* ISTATE state vector containing the library parameters.
* MATNO reaction material indices.
*
*Parameters: input/output
* KPAX complete reaction type matrix.
* BPAX complete branching ratio matrix.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER MAXR,NEL,NBESP,NSTATE,ITNAM(3,NEL),ISTATE(NSTATE),
> MATNO(NEL),KPAX(NEL+MAXR,NEL)
CHARACTER NMDEPL(MAXR)*8
REAL BPAX(NBESP,NEL+MAXR,NEL)
*----
* LOCAL VARIABLES
*----
INTEGER NDFP,JSO,ISO,ITR,NUNTIL,NHEAVY,NDFI,IUNTIL,NEW,
> NLIGHT,NOTHER,NIFP,KPAX0,KPAX1,KPAX2,IEL,JEL,
> KREAC,II,NPAR,ICOUNT,NREAC,NSTABL
REAL SUMFI
CHARACTER TEXT12*12
LOGICAL LOGPF
*----
* INTERNAL PARAMETERS
*----
INTEGER KDECAY,KFISSP,KCAPTU,KFISSI,KHEAT
PARAMETER (KDECAY=1,KFISSP=2,KCAPTU=3,KFISSI=6,KHEAT=9)
*----
* COMPUTE NUMBER OF DIRECT FISSION PRODUCT
* NUMBER DIRECT FISSION PRODUCT
*
* SET MATNO TO LIGHT (-KFISSP) FOR FISSION PRODUCT
* SET MATNO TO HEAVY (-KFISSI) FOR FISSILE ISOTOPES
* SET MATNO TO STABLE (-KHEAT) FOR STABLE ISOTOPES PRODUCING ENERGY
* SET MATNO TO 0 TO UNUSED ISOTOPES
*----
NDFI=0
NDFP=0
DO 100 JSO=NEL,1,-1
MATNO(JSO)=0
KPAX0=KPAX(NEL+KDECAY,JSO)
KPAX1=KPAX(NEL+KFISSP,JSO)
KPAX2=KPAX(NEL+KCAPTU,JSO)
IF((KPAX0.EQ.-9999).OR.(KPAX1.EQ.-9999).OR.(KPAX2.EQ.-9999))
> THEN
* JSO IS A STABLE ISOTOPE
MATNO(JSO)=-KHEAT
DO 222 ITR=1,MAXR
IF(BPAX(1,NEL+ITR,JSO).GT.0.0) THEN
KPAX(NEL+ITR,JSO)=2
KPAX(JSO,:NEL)=0
ENDIF
222 CONTINUE
GO TO 100
ELSE IF(KPAX1.LT.0) THEN
* JSO IS A FISSION PRODUCT
MATNO(JSO)=-KFISSP
NIFP=0
DO 101 ISO=NEL,1,-1
IF((KPAX(JSO,ISO).EQ.KFISSP).AND.
> (BPAX(1,JSO,ISO).GT.0.0)) THEN
NIFP=NIFP+1
ENDIF
101 CONTINUE
IF(NIFP.GT.0) THEN
NDFP=NDFP+1
KPAX(NEL+KFISSP,JSO)=5+100*NDFP
ELSE
KPAX(NEL+KFISSP,JSO)=5
ENDIF
ELSE IF(KPAX1.GT.0) THEN
* JSO IS A FISSILE ISOTOPE
MATNO(JSO)=-KFISSI
SUMFI=0.0
DO 221 ISO=1,NEL
IF(KPAX(ISO,JSO).EQ.KFISSP) SUMFI=SUMFI+BPAX(1,ISO,JSO)
221 CONTINUE
IF(SUMFI.GT.0.0) THEN
NDFI=NDFI+1
KPAX(NEL+KFISSP,JSO)=4+100*NDFI
ELSE
KPAX(NEL+KFISSP,JSO)=3
ENDIF
ENDIF
100 CONTINUE
*----
* CHECK IF THE DEPLETION CHAIN IS COHERENT
*----
DO 200 IEL=1,NEL
DO 210 JEL=1,NEL
KREAC=KPAX(JEL,IEL)
IF(KREAC.EQ.KFISSP) THEN
IF(MOD(KPAX(NEL+KFISSP,JEL),100).NE.5) THEN
WRITE(TEXT12,'(3A4)') (ITNAM(II,JEL),II=1,3)
CALL XABORT('LIBWET: SON '//TEXT12//' IS NOT A FISSION '//
> 'PRODUCT')
ENDIF
IF((KPAX(NEL+KFISSP,IEL).NE.3).AND.
> (MOD(KPAX(NEL+KFISSP,IEL),100).NE.4)) THEN
WRITE(TEXT12,'(3A4)') (ITNAM(II,IEL),II=1,3)
CALL XABORT('LIBWET: PARENT '//TEXT12//' IS NOT A FISSI'//
> 'LE ISOTOPE')
ENDIF
ELSE IF(KREAC.NE.0) THEN
IF(KPAX(NEL+KREAC,IEL).EQ.0) THEN
WRITE(TEXT12,'(3A4)') (ITNAM(II,IEL),II=1,3)
CALL XABORT('LIBWET: PARENT '//TEXT12//' DOES NOT DEPL'//
> 'ETE VIA REACTION '//NMDEPL(KREAC))
ENDIF
ENDIF
210 CONTINUE
200 CONTINUE
*----
* SET MATNO TO OTHER (-KDECAY) FOR ISOTOPES WITH DAUGHTER OR FATHER
* AND FOR ISOTOPES WITH DECAY
*----
DO 112 ISO=1,NEL
IF(MATNO(ISO).EQ.0) THEN
IF(KPAX(NEL+KDECAY,ISO).GT.0 .OR.
> KPAX(NEL+KCAPTU,ISO).GT.0) THEN
MATNO(ISO)=-KDECAY
GO TO 115
ENDIF
DO 113 JSO=1,NEL
IF(KPAX(JSO,ISO).GT.0) THEN
MATNO(ISO)=-KDECAY
GO TO 115
ELSE IF(KPAX(ISO,JSO).GT.0) THEN
MATNO(ISO)=-KDECAY
GO TO 115
ENDIF
113 CONTINUE
ENDIF
115 CONTINUE
112 CONTINUE
*----
* SET MATNO TO STABLE (-KHEAT) FOR OTHER ISOTOPES PRODUCING ENERGY
*----
DO 212 ISO=1,NEL
IF(MATNO(ISO).EQ.0) THEN
DO 213 ITR=2,MAXR
IF(BPAX(1,NEL+ITR,ISO).NE.0.0) THEN
MATNO(ISO)=-KHEAT
GO TO 215
ENDIF
213 CONTINUE
ENDIF
215 CONTINUE
212 CONTINUE
*----
* COMPUTE THE MAXIMUM NUMBER OF PARENTS FOR AN ISOTOPE
*----
NPAR=0
DO 116 ISO=1,NEL
ICOUNT=0
DO 114 JSO=1,NEL
LOGPF=((MATNO(ISO).EQ.-KFISSP).AND.(MATNO(JSO).EQ.-KFISSI))
IF((BPAX(1,ISO,JSO).NE.0.0).AND.(.NOT.LOGPF)) ICOUNT=ICOUNT+1
114 CONTINUE
NPAR=MAX(NPAR,ICOUNT)
116 CONTINUE
*----
* COMPUTE THE MAXIMUM NUMBER OF DEPLETION REACTIONS
*----
NREAC=4
DO 118 ISO=1,NEL
DO 117 ITR=1,MAXR
IF(KPAX(NEL+ITR,ISO).NE.0) NREAC=MAX(NREAC,ITR)
117 CONTINUE
118 CONTINUE
*----
* SET MATNO TO HEAVY (-KFISSI) FOR ISOTOPES RESULTING FROM DECAY
* OR CAPTURE OF HEAVY ISOTOPE UNTIL ALL HEAVY ISOTOPES IDENTIFIED
*----
NUNTIL=NEL
NHEAVY=0
DO 120 IUNTIL=1,NUNTIL
NEW=0
DO 121 ISO=NEL,1,-1
IF(MATNO(ISO).EQ.-KFISSI) THEN
NHEAVY=NHEAVY+1
DO 122 JSO=NEL,1,-1
IF(MATNO(JSO).NE.-KFISSI) THEN
IF((KPAX(JSO,ISO).NE.0).AND.(KPAX(JSO,ISO).NE.KFISSP))
> THEN
NEW=NEW+1
MATNO(JSO)=-KFISSI
ENDIF
ENDIF
122 CONTINUE
ENDIF
121 CONTINUE
IF(NEW.EQ.0) GO TO 125
NHEAVY=0
120 CONTINUE
125 CONTINUE
*----
* SET MATNO TO LIGHT (-KFISSP) FOR ISOTOPES RESULTING FROM DECAY
* OR CAPTURE OF LIGHT ISOTOPE UNTIL ALL LIGHT ISOTOPES IDENTIFIED
*----
NUNTIL=NUNTIL-NHEAVY
NLIGHT=0
DO 130 IUNTIL=1,NUNTIL
NEW=0
DO 131 ISO=NEL,1,-1
IF(MATNO(ISO).EQ.-KFISSP) THEN
NLIGHT=NLIGHT+1
DO 132 JSO=NEL,1,-1
IF(MATNO(JSO).NE.-KFISSI.AND.
> MATNO(JSO).NE.-KFISSP) THEN
IF(KPAX(JSO,ISO).NE.0) THEN
NEW=NEW+1
MATNO(JSO)=-KFISSP
ENDIF
ENDIF
132 CONTINUE
ENDIF
131 CONTINUE
IF(NEW.EQ.0) GO TO 135
NLIGHT=0
130 CONTINUE
135 CONTINUE
*----
* SET MATNO TO OTHER (-KDECAY) FOR ISOTOPES RESULTING FROM DECAY
* OR CAPTURE OF OTHER ISOTOPE UNTIL ALL OTHER ISOTOPES IDENTIFIED
*----
NUNTIL=NUNTIL-NLIGHT
NOTHER=0
DO 140 IUNTIL=1,NUNTIL
NEW=0
DO 141 ISO=NEL,1,-1
IF(MATNO(ISO).EQ.-KDECAY) THEN
NOTHER=NOTHER+1
DO 142 JSO=NEL,1,-1
IF(MATNO(JSO).NE.-KFISSI.AND.
> MATNO(JSO).NE.-KFISSP.AND.
> MATNO(JSO).NE.-KDECAY) THEN
IF(KPAX(JSO,ISO).NE.0) THEN
NEW=NEW+1
MATNO(JSO)=-KDECAY
ENDIF
ENDIF
142 CONTINUE
ENDIF
141 CONTINUE
IF(NEW.EQ.0) GO TO 145
NOTHER=0
140 CONTINUE
145 CONTINUE
*----
* COUNT THE NUMBER OF STABLE ISOTOPES AND SET TO ZERO THEIR RADIOACTIVE
* DECAY CONSTANT.
*----
NSTABL=0
DO 150 ISO=1,NEL
IF(MATNO(ISO).EQ.-KHEAT) THEN
NSTABL=NSTABL+1
BPAX(:,NEL+KDECAY,ISO)=0.0
ENDIF
150 CONTINUE
*
ISTATE(:NSTATE)=0
ISTATE(1)=NHEAVY+NLIGHT+NOTHER+NSTABL
ISTATE(2)=NDFI
ISTATE(3)=NDFP
ISTATE(4)=NHEAVY
ISTATE(5)=NLIGHT
ISTATE(6)=NOTHER
ISTATE(7)=NSTABL
ISTATE(8)=NREAC
ISTATE(9)=NPAR
ISTATE(10)=NBESP
*----
* RETURN
*----
RETURN
END
|
