summaryrefslogtreecommitdiff
path: root/Dragon/src/LIBEPR.f
blob: ad893cc0765063e5bf98a9c761d47514490678ac (plain)
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
*DECK LIBEPR
      SUBROUTINE LIBEPR(IMPX,NBESP,NDEPL,NSTABL,NDFI,NDFP,NREAC,NPAR,
     >                  HNADPL,NMDEPL,IDR,RER,RRD,KPAR,BPAR,YIELD,IZAE)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Print isotopic depletion chain.
* 
*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
* IMPX    print parameter:
*         .ge.5 nice formatted printout;
*         .ge.10 dragon input formatted;
*         .ge.20 dragon/APOLIB-2 input formatted.
* NBESP   number of energy-dependent fission yield matrices.
* NDEPL   number of depleting isotopes.
* NSTABL  number of non-depleting isotopes producing energy.
* NDFI    number of direct fissile isotopes.
* NDFP    number of direct fission product.
* NREAC   maximum number of depletion reaction in the depletion chain.
* NPAR    maximum number of parent isopopes from decay and 
*         neutron-induced reactions.
* HNADPL  reactive isotope names in chain.
* NMDEPL  names of used depletion reactions
*         (NMDEPL(1)='DECAY'; NMDEPL(2)='NFTOT';
*          NMDEPL(3)='NG'   ; NMDEPL(4)='N2N';
*          etc.).
* IDR     DEPLETE-REAC matrix (reaction identifiers).
* RER     DEPLETE-ENER matrix (MeV/reaction values).
* RRD     DEPLETE-DECA vector (decay constant values).
* KPAR    PRODUCE-REAC matrix (production identifiers).
* BPAR    PRODUCE-RATE matrix (branching ratios).
* YIELD   fission product yield matrix.
* IZAE    6-digit nuclide identifier:
*         atomic number z*10000 (digits) + mass number a*10 +
*         energy state (0 = ground state, 1 = first state, etc.).
*
*-----------------------------------------------------------------------
*
      IMPLICIT NONE
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER    IMPX,NBESP,NDEPL,NSTABL,NDFI,NDFP,NREAC,NPAR,
     >           HNADPL(3,NDEPL),NMDEPL(2,NREAC),IDR(NREAC,NDEPL),
     >           KPAR(NPAR,NDEPL),IZAE(NDEPL)
      REAL       RER(NREAC,NDEPL),RRD(NDEPL),BPAR(NPAR,NDEPL),
     >           YIELD(NBESP,NDFI,NDFP)
*----
*  LOCAL VARIABLES
*----
      INTEGER    I,J,ISO,KRD,KRF,ISOF,IPF,KRC,IOF,KROTH,NBPROD,IPAR,
     >           JSO,KREAC,ISOFP,IOFS,IREAC,NBFIS,IOUT,IBESP
      REAL       RATD1,RATD2,RATF,RATC,RREAC
      REAL, ALLOCATABLE, DIMENSION(:) :: RATFS
*----
*  INTERNAL PARAMETERS
*----
      REAL        CNVDAY
      PARAMETER  (IOUT=6,CNVDAY=1.0E+8/86400.0)
      CHARACTER   NONOFF(2)*3
      CHARACTER   LINE*130
      SAVE        NONOFF
      DATA        NONOFF/'NO ','YES'/
*----
*  PRINT THE DEPLETION CHAIN.
*----
      IF(IMPX.GE.5) THEN
        WRITE(IOUT,6000) ((NMDEPL(J,I),J=1,2),I=4,NREAC)
        WRITE(IOUT,6010)
        ALLOCATE(RATFS(NBESP))
        DO 40 ISO=1,NDEPL
          KRD=1
          RATD1=0.0
          RATD2=0.0
          IF(IDR(1,ISO).NE.0) THEN
            RATD1=CNVDAY/RRD(ISO)
            RATD2=RER(1,ISO)
            KRD=2
          ENDIF
          KRF=1
          RATF=0.0
          RATFS(:)=0.0
          IF(MOD(IDR(2,ISO),100).EQ.3) THEN
            KRF=2
            RATF=RER(2,ISO)
          ELSE IF(MOD(IDR(2,ISO),100).EQ.4) THEN
            KRF=2
            RATF=RER(2,ISO)
            ISOF=IDR(2,ISO)/100
            IF(ISOF.NE.0) THEN
              DO 20 IBESP=1,NBESP
                DO 10 IPF=1,NDFP
                  RATFS(IBESP)=RATFS(IBESP)+YIELD(IBESP,ISOF,IPF)
  10            CONTINUE
  20          CONTINUE
            ENDIF
          ENDIF
          KRC=1
          RATC=0.0
          IF(IDR(3,ISO).NE.0) THEN
            RATC=RER(3,ISO)
            KRC=2
          ENDIF
          LINE=' '
*----
*  WRITE ISOTOPE PROPERTIES
*----
          WRITE(LINE(:18),'(I5,1X,2A4,1X,A3)') ISO,HNADPL(1,ISO),
     >    HNADPL(2,ISO),NONOFF(KRD)
          IF(KRD.EQ.2) WRITE(LINE(19:42),'(1P,2E12.4)') RATD1,RATD2
          WRITE(LINE(45:47),'(A3)') NONOFF(KRF)
          IF(KRF.EQ.2) WRITE(LINE(48:71),'(1P,2E12.4)') RATF,RATFS(1)
          WRITE(LINE(74:76),'(A3)') NONOFF(KRC)
          IF(KRC.EQ.2) WRITE(LINE(77:88),'(1P,E12.4)') RATC
          IOF=91
          DO 30 I=4,NREAC
            KROTH=1
            IF(IDR(I,ISO).GT.0) KROTH=2
            IF(IOF+7.GT.130) THEN
              WRITE(IOUT,'(1X,A)') LINE
              IOF=91
              LINE=' '
            ENDIF
            WRITE(LINE(IOF:IOF+7),'(A3,5X)') NONOFF(KROTH)
            IOF=IOF+8
  30      CONTINUE
          WRITE(IOUT,'(1X,A)') LINE
  40    CONTINUE
        DEALLOCATE(RATFS)
*----
*  WRITE PARENTS FROM ALL REACTION EXCEPT FISSION
*----
        WRITE(IOUT,7000)
        DO 70 ISO=1,NDEPL-NSTABL
          NBPROD=0
          DO 50 IPAR=1,NPAR
            JSO=KPAR(IPAR,ISO)/100
            KREAC=MOD(KPAR(IPAR,ISO),100)
            RREAC=BPAR(IPAR,ISO)
            IF(JSO.GT.0) THEN
              IF((KREAC.LE.0).OR.(KREAC.GT.NREAC)) CALL XABORT('LIBEPR'
     >        //': INDALID REACTION INDEX')
              NBPROD=NBPROD+1
              IF(NBPROD.EQ.1) THEN
                WRITE(IOUT,6012) ISO,HNADPL(1,ISO),HNADPL(2,ISO),
     >          NMDEPL(1,KREAC),NMDEPL(2,KREAC),JSO,HNADPL(1,JSO),
     >          HNADPL(2,JSO),RREAC
              ELSE
                WRITE(IOUT,6013) NMDEPL(1,KREAC),NMDEPL(2,KREAC),
     >          JSO,HNADPL(1,JSO),HNADPL(2,JSO),RREAC
              ENDIF
            ENDIF
  50      CONTINUE
*----
*  WRITE PARENTS FROM FISSION IF REQUIRED
*----
          IF(MOD(IDR(2,ISO),100).EQ.2) THEN
            GO TO 70
          ELSE IF(MOD(IDR(2,ISO),100).EQ.5) THEN
            ISOFP=IDR(2,ISO)/100
            IF(ISOFP.GT.NDFP)
     >        CALL XABORT('LIBEPR: INVALID FISSION PRODUCT NUMBER')
            DO 60 JSO=1,NDEPL
              IF(MOD(IDR(2,JSO),100).EQ.4) THEN
                ISOF=IDR(2,JSO)/100
                IF(ISOF.GT.NDFI) THEN
                  CALL XABORT('LIBEPR: INVALID FISSILE NUMBER')
                ELSE IF(ISOF.GT.0) THEN
                  RREAC=YIELD(1,ISOF,ISOFP)
                  IF(RREAC.GT.0.0) THEN
                    NBPROD=NBPROD+1
                    IF(NBPROD.EQ.1) THEN
                      WRITE(IOUT,6012) ISO,HNADPL(1,ISO),HNADPL(2,ISO),
     >                NMDEPL(1,2),NMDEPL(2,2),JSO,HNADPL(1,JSO),
     >                HNADPL(2,JSO),RREAC
                    ELSE
                      WRITE(IOUT,6013) NMDEPL(1,2),NMDEPL(2,2),JSO,
     >                HNADPL(1,JSO),HNADPL(2,JSO),RREAC
                    ENDIF
                  ENDIF
                ENDIF
              ENDIF
  60        CONTINUE
          ENDIF
  70    CONTINUE
      ENDIF
*
      IF(IMPX.GE.10) THEN
        WRITE(IOUT,'(/1X,4HDEPL,I6,6H CHAIN)') NDEPL
        DO 330 ISO=1,NDEPL
        LINE=' '
        WRITE(LINE(:17),'(1H'',2A4,1H'',I7)') HNADPL(1,ISO),
     >  HNADPL(2,ISO),IZAE(ISO)
        IOFS=18
        IF(IDR(1,ISO).NE.0) THEN
          WRITE(LINE(IOFS:IOFS+20),'(1X,5HDECAY,1P,E14.6)') RRD(ISO)
          IOFS=IOFS+20
*          IF(RER(1,ISO).NE.0.0) THEN
*            WRITE(LINE(IOFS:IOFS+13),'(1X,1P,E12.5)') RER(1,ISO)
*            IOFS=IOFS+13
*          ENDIF
        ENDIF
        IF((MOD(IDR(2,ISO),100).EQ.3).OR.(MOD(IDR(2,ISO),100).EQ.4))
     >  THEN
          WRITE(LINE(IOFS:IOFS+15),'(1X,5HNFTOT,F9.4)') RER(2,ISO)
          IOFS=IOFS+15
        ENDIF
        IF(IDR(3,ISO).NE.0) THEN
          WRITE(LINE(IOFS:IOFS+17),'(1X,2HNG,1P,E14.6)') RER(3,ISO)
          IOFS=IOFS+17
        ENDIF
        DO 80 IREAC=4,NREAC
        IF(IDR(IREAC,ISO).NE.0) THEN
          IF(IOFS+9.GT.71) THEN
             WRITE(IOUT,'(1X,A)') LINE
             LINE=' '
             IOFS=18
          ENDIF
          WRITE(LINE(IOFS:IOFS+9),'(1X,2A4)') NMDEPL(1,IREAC),
     >    NMDEPL(2,IREAC)
          IOFS=IOFS+9
        ENDIF
  80    CONTINUE
        IF(ISO.GT.NDEPL-NSTABL) THEN
          WRITE(LINE(IOFS:IOFS+7),'(7H STABLE)')
          IOFS=IOFS+7
          IF(IOFS.GT.71) CALL XABORT('LIBEPR: LINE OVERFLOW.')
        ENDIF
        IF(LINE.NE.' ') WRITE(IOUT,'(1X,A)') LINE
        LINE=' '
        IOFS=7
        NBPROD=0
        DO 310 IPAR=1,NPAR
          JSO=KPAR(IPAR,ISO)/100
          KREAC=MOD(KPAR(IPAR,ISO),100)
          RREAC=BPAR(IPAR,ISO)
          IF(JSO.GT.0) THEN
            NBPROD=NBPROD+1
            IF(NBPROD.EQ.1) WRITE(LINE(2:6),'(5HFROM )')
            IF(IOFS+34.GE.71) THEN
               WRITE(IOUT,'(1X,A)') LINE
               LINE=' '
               IOFS=7
            ENDIF
            WRITE(LINE(IOFS:IOFS+34),6100)
     >      NMDEPL(1,KREAC),NMDEPL(2,KREAC),RREAC,HNADPL(1,JSO),
     >      HNADPL(2,JSO)
            IOFS=IOFS+34
          ENDIF
 310    CONTINUE
        IF(LINE.NE.' ') WRITE(IOUT,'(1X,A)') LINE
        LINE=' '
        IOFS=16
        NBFIS=0
        IF(MOD(IDR(2,ISO),100).EQ.2) THEN
          GO TO 330
        ELSE IF(MOD(IDR(2,ISO),100).EQ.5) THEN
          ISOFP=IDR(2,ISO)/100
          DO 320 JSO=1,NDEPL
            IF(MOD(IDR(2,JSO),100).EQ.4) THEN
              ISOF=IDR(2,JSO)/100
              IF(ISOF.GT.0) THEN
                RREAC=YIELD(1,ISOF,ISOFP)
                IF(RREAC.GT.0.0) THEN
                  NBPROD=NBPROD+1
                  IF(NBPROD.EQ.1) WRITE(LINE(2:6),'(5HFROM )')
                  NBFIS=NBFIS+1
                  IF(NBFIS.EQ.1) WRITE(LINE(7:14),'(A8)') 'NFTOT   '
                  WRITE(LINE(IOFS:IOFS+24),6101)
     >            RREAC,HNADPL(1,JSO),HNADPL(2,JSO)
                  IOFS=IOFS+24
                  IF(IOFS.GE.60) THEN
                     WRITE(IOUT,'(1X,A)') LINE
                     LINE=' '
                     IOFS=16
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
 320      CONTINUE
          IF(LINE.NE.' ') WRITE(IOUT,'(1X,A)') LINE
        ENDIF
 330    CONTINUE
        WRITE(IOUT,'(9H ENDCHAIN)')
      ENDIF
*
      IF(IMPX.GE.500) THEN
        WRITE(IOUT,'(/1X,33HDEPL LIB: APLIB2 FIL: XXXXX CHAIN)')
        DO 350 ISO=1,NDEPL-NSTABL
        LINE=' '
        WRITE(LINE(:8),'(2A4)') HNADPL(1,ISO),HNADPL(2,ISO)
        IOFS=14
        NBPROD=0
        DO 340 IPAR=1,NPAR
          JSO=KPAR(IPAR,ISO)/100
          KREAC=MOD(KPAR(IPAR,ISO),100)
          RREAC=BPAR(IPAR,ISO)
          IF(JSO.GT.0) THEN
            NBPROD=NBPROD+1
            IF(NBPROD.EQ.1) WRITE(LINE(10:13),'(4HFROM)')
            WRITE(LINE(IOFS:IOFS+29),'(1X,2A4,1P,E11.4,1X,2A4)')
     >      NMDEPL(1,KREAC),NMDEPL(2,KREAC),RREAC,HNADPL(1,JSO),
     >      HNADPL(2,JSO)
            IOFS=IOFS+29
            IF(IOFS.GE.71) THEN
               WRITE(IOUT,'(1X,A)') LINE
               LINE=' '
               IOFS=14
            ENDIF
          ENDIF
 340    CONTINUE
        IF(LINE.NE.' ') WRITE(IOUT,'(1X,A)') LINE
 350    CONTINUE
        WRITE(IOUT,'(9H ENDCHAIN)')
      ENDIF
*----
*  RETURN FROM LIBEPR
*----
      RETURN
*----
*  FORMAT
*----
 6000 FORMAT(/' DEPLETION CHAIN DATA'/' --------------------'//
     >       43X,'DEPLETION REACTIONS'/
     >       '                -----------------------------',
     >       '--------------------------------------------'/
     >       '   ISOTOPE      ......DECAY................  ..........',
     >       '.FISSION.........  ......NG.......  ',10A4:/(91X,10A4))
 6010 FORMAT('   NB. NAME         MLIFE(DAYS) ENERGY(MEV)      ENERGY',
     >       '(MEV) TOTAL YIELD      ENERGY(MEV)')
 6012 FORMAT(I5,1X,2A4,4X,2A4,3X,I5,1X,2A4,1X,1P,E14.6)
 6013 FORMAT(18X,2A4,3X,I5,1X,2A4,1X,1P,E14.6)
 6100 FORMAT(2A4,1X,1P,E13.6,2H ',2A4,2H' )
 6101 FORMAT(1P,E13.6,2H ',2A4,1H')
 7000 FORMAT(//27X,'PRODUCTION REACTIONS'/
     >       18X,'---------------------------------------'/
     >       '  ISOTOPE         REACTION     ISOTOPE',14X,'YIELD'/
     >       '  NB. NAME                     NB. NAME    ')
      END