summaryrefslogtreecommitdiff
path: root/Donjon/src/T16ENE.f
blob: b505bf135ed6bfcfce68d97df8bddc50319e0c95 (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
*DECK T16ENE
      SUBROUTINE T16ENE(IPRINT,MXGRP ,NG    ,NGCOND,NGMTR ,NGREAC,
     >                  NGCCPO,IFGCND,IFGMTR,IFGEDI,ENECPO,ENET16,
     >                  VELMTR)
*
*----
*
*Purpose:
*  Generate and analyse energy structure.
*
*Author(s): 
* G. Marleau
*
*Parameters: input
* IPRINT  print level.
* MXGRP   maximum number of groups.
* NG      number of groups in library.
* NGCOND  number of condensed groups.
* NGMTR   numbre of main transport groups.
* NGREAC  numbre of edit groups.
*
*Parameters: input/output
* NGCCPO  numbre of edit groups.
* IFGCND  reference/exit condensation few groups.
* IFGMTR  reference/exit main transport few groups.
* IFGEDI  reference/exit edit few groups.
* ENECPO  final energy group structure for CPO.
*
*Parameters: output
* ENET16  energy group structure for tape16.
* VELMTR  velocity for main transport.
*
*----
*
      IMPLICIT         NONE
      INTEGER          IPRINT,MXGRP,NG,NGCOND,NGMTR,NGREAC,NGCCPO
      INTEGER          IFGCND(MXGRP),IFGMTR(MXGRP),
     >                 IFGEDI(MXGRP)
      REAL             ENECPO(MXGRP+1),ENET16(MXGRP+1),
     >                 VELMTR(MXGRP)
*----
*  LOCAL VARIABLES
*  FOR AVERAGED NEUTRON VELOCITY
*  V=SQRT(2*ENER/M)=SQRT(2/M)*SQRT(ENER)
*  SQFMAS=SQRT(2/M) IN CM/S/SQRT(EV) FOR V IN CM/S AND E IN EV
*        =SQRT(2*1.602189E-19(J/EV)* 1.0E4(CM2/M2) /1.67495E-27 (KG))
*        =1383155.30602 CM/S/SQRT(EV)
*----
      INTEGER          IOUT,MGELIB,MGWLIB
      CHARACTER        NAMSBR*6
      REAL             SQFMAS,PRECIS
      PARAMETER       (IOUT=6,MGELIB=89,MGWLIB=69,
     >                 NAMSBR='T16ENE',SQFMAS=1383155.30602,
     >                 PRECIS=1.0E-5)
      INTEGER          IGR,IGC,IGD,IGF
      REAL             EAVG
      INTEGER, ALLOCATABLE, DIMENSION(:) :: IFGCPO
      REAL, ALLOCATABLE, DIMENSION(:) :: VELEDI,VELT16
*----
*  DATA
*----
      REAL             ENEELB(MGELIB+1),ENEWLB(MGWLIB+1)
      SAVE             ENEELB,ENEWLB
      DATA             ENEELB
     >/1.0000E+7,7.7880E+6,6.0653E+6,4.7237E+6,3.6788E+6,2.8650E+6,
     > 2.2313E+6,1.7377E+6,1.3534E+6,1.0540E+6,8.2085E+5,6.3928E+5,
     > 4.9787E+5,3.8774E+5,3.0197E+5,2.3518E+5,1.8316E+5,1.4264E+5,
     > 1.1109E+5,8.6517E+4,6.7379E+4,4.0868E+4,2.4788E+4,1.5034E+4,
     > 9.1188E+3,5.5308E+3,3.3546E+3,2.0347E+3,1.2341E+3,7.4852E+2,
     > 4.5400E+2,2.7536E+2,1.6702E+2,1.3007E+2,1.0130E+2,7.8893E+1,
     > 6.1442E+1,4.7851E+1,3.7267E+1,2.9023E+1,2.2603E+1,1.7603E+1,
     > 1.3710E+1,1.0677E+1,8.3153E+0,6.4760E+0,5.0435E+0,4.0000E+0,
     > 3.3000E+0,2.6000E+0,2.1000E+0,1.5000E+0,1.3000E+0,1.1500E+0,
     > 1.1230E+0,1.0970E+0,1.0710E+0,1.0450E+0,1.0200E+0,9.9600E-1,
     > 9.7200E-1,9.5000E-1,9.1000E-1,8.5000E-1,7.8000E-1,6.2500E-1,
     > 5.0000E-1,4.0000E-1,3.5000E-1,3.2000E-1,3.0000E-1,2.8000E-1,
     > 2.5000E-1,2.2000E-1,1.8000E-1,1.4000E-1,1.0000E-1,8.0000E-2,
     > 6.7000E-2,5.8000E-2,5.0000E-2,4.2000E-2,3.5000E-2,3.0000E-2,
     > 2.5000E-2,2.0000E-2,1.5000E-2,1.0000E-2,5.0000E-3,2.0000E-4/
      DATA             ENEWLB
     >/1.0000E+7,6.0655E+6,3.6790E+6,2.2310E+6,1.3530E+6,8.2100E+5,
     > 5.0000E+5,3.0250E+5,1.8300E+5,1.1100E+5,6.7340E+4,4.0850E+4,
     > 2.4780E+4,1.5030E+4,9.1180E+3,5.5300E+3,3.5191E+3,2.2394E+3,
     > 1.4251E+3,9.0690E+2,3.6726E+2,1.4873E+2,7.5501E+1,4.8052E+1,
     > 2.7700E+1,1.5968E+1,9.8770E+0,4.0000E+0,3.3000E+0,2.6000E+0,
     > 2.1000E+0,1.5000E+0,1.3000E+0,1.1500E+0,1.1230E+0,1.0970E+0,
     > 1.0710E+0,1.0450E+0,1.0200E+0,9.9600E-1,9.7200E-1,9.5000E-1,
     > 9.1000E-1,8.5000E-1,7.8000E-1,6.2500E-1,5.0000E-1,4.0000E-1,
     > 3.5000E-1,3.2000E-1,3.0000E-1,2.8000E-1,2.5000E-1,2.2000E-1,
     > 1.8000E-1,1.4000E-1,1.0000E-1,8.0000E-2,6.7000E-2,5.8000E-2,
     > 5.0000E-2,4.2000E-2,3.5000E-2,3.0000E-2,2.5000E-2,2.0000E-2,
     > 1.5000E-2,1.0000E-2,5.0000E-3,1.0000E-5/
*----
*  STORE ORIGINAL GROUP STRUCTURE IN ENET16
*----
      ALLOCATE(IFGCPO(MXGRP),VELEDI(MXGRP),VELT16(MXGRP))
      IF(NG .EQ. MGELIB) THEN
        ENET16(1)=ENEELB(1)
        DO IGR=2,MGELIB+1
          ENET16(IGR)=ENEELB(IGR)
          EAVG=SQRT(ENET16(IGR)*ENET16(IGR-1))
          VELT16(IGR-1)=SQFMAS*SQRT(EAVG)
        ENDDO
      ELSE IF(NG .EQ. MGWLIB) THEN
        ENET16(1)=ENEWLB(1)
        DO IGR=2,MGWLIB+1
          ENET16(IGR)=ENEWLB(IGR)
          EAVG=SQRT(ENET16(IGR)*ENET16(IGR-1))
          VELT16(IGR-1)=SQFMAS*SQRT(EAVG)
        ENDDO
      ELSE
        CALL XABORT(NAMSBR//
     >  ': INVALID TAPE16 GROUP STRUCTURE')
      ENDIF
      IF(IPRINT .GE. 10) THEN
        WRITE(IOUT,6000) NAMSBR
        WRITE(IOUT,6010) NG,NGMTR,NGREAC,NGCOND,NGCCPO
        WRITE(IOUT,6030) (ENET16(IGC),IGC=1,NG+1)
        WRITE(IOUT,6040) (VELT16(IGC),IGC=1,NG)
      ENDIF
*----
*  COMPUTE AVERAGED NEUTRON GROUP VELOCITY
*  AVERAGED NEUTRON ENERGY ENER=SQRT(E(G+1)*E(G))
*  V=SQRT(2*ENER/M)=SQRT(2/M)*SQRT(ENER)
*   =SQFMAS*SQRT(ENER)
*----
      IF(NGMTR .GT. 0) THEN
        IGF=1
        DO IGR=1,NGMTR
          IGD=IGF
          IGF=IFGMTR(IGR)+1
          EAVG=SQRT(ENET16(IGD)*ENET16(IGF))
          VELMTR(IGR)=SQFMAS*SQRT(EAVG)
        ENDDO
      ENDIF
      IF(IPRINT .GE. 10 .AND. NGMTR .GT. 0) THEN
        WRITE(IOUT,6020) (IFGMTR(IGC),IGC=1,NGMTR)
        WRITE(IOUT,6031) ENET16(1),
     >  (ENET16(IFGMTR(IGC)+1),IGC=1,NGMTR)
        WRITE(IOUT,6041) (VELMTR(IGC),IGC=1,NGMTR)
      ENDIF
      IF(NGREAC .GT. 0) THEN
        IGF=1
        DO IGR=1,NGREAC
          IGD=IGF
          IGF=IFGEDI(IGR)+1
          EAVG=SQRT(ENET16(IGD)*ENET16(IGF))
          VELEDI(IGR)=SQFMAS*SQRT(EAVG)
        ENDDO
*----
*    TEST IF CONDENSATION STRUCTURE PROVIDED BY IFGEDI
*    COMPATIBLE WITH IFGMTR
*----
        IF(NGMTR .GT. 0) THEN
          DO IGC=1,NGREAC
            DO IGR=IGC,NGMTR
              IF(IFGEDI(IGC) .EQ. IFGMTR(IGR)) THEN
                GO TO 105
              ENDIF
            ENDDO
            CALL XABORT(NAMSBR//
     >      ': IFGEDI AND IFGMTR NOT COMPATIBLE')
 105        CONTINUE
          ENDDO
        ENDIF
      ENDIF
      IF(IPRINT .GE. 10 .AND. NGREAC .GT. 0) THEN
        WRITE(IOUT,6021) (IFGEDI(IGC),IGC=1,NGREAC)
        WRITE(IOUT,6032) ENET16(1),
     >  (ENET16(IFGEDI(IGC)+1),IGC=1,NGREAC)
        WRITE(IOUT,6042) (VELEDI(IGC),IGC=1,NGREAC)
      ENDIF
*----
*  IF NGCCPO > 0 FIND IFGCPO FROM ENECPO
*----
      IF(NGCCPO .GT. 0) THEN
        IF(ABS(ENECPO(1)-ENET16(1)) .GT. PRECIS ) CALL XABORT(NAMSBR//
     >  ': ENECPO(1) SHOULD BE IDENTICAL TO ENET16(1)')
        DO IGC=2,NGCCPO+1
          DO IGR=IGC,NG+1
            IF(ABS(ENECPO(IGC)-ENET16(IGR)) .LT. PRECIS ) THEN
              IFGCPO(IGC-1)=IGR-1
              GO TO 115
            ENDIF
          ENDDO
 115      CONTINUE
        ENDDO
        IF(NGCOND .GT. 0) THEN
*----
*  IF NGCOND > 0
*  IFGCPO AND IFGCND NUST BE IDENTICAL
*----
          IF(NGCCPO .NE. NGCOND) CALL XABORT(NAMSBR//
     >    ': NGCCPO AND NGCOND MUST BE IDENTICAL')
          DO IGC=1,NGCCPO
            IF(IFGCPO(IGC) .NE. IFGCND(IGC))
     >      CALL XABORT(NAMSBR//
     >      ': IFGCPO AND IFGCND MUST BE IDENTICAL')
          ENDDO
        ENDIF
      ELSE
*----
*  IF NGCCPO =0
*----
        IF(NGCOND .GT. 0) THEN
*----
*  IF NGCOND > 0
*  IFGCPO = IFGCND
*----
          NGCCPO=NGCOND
          ENECPO(1)=ENET16(1)
          DO IGC=1,NGCCPO
            IFGCPO(IGC)=IFGCND(IGC)
            ENECPO(IGC+1)=ENET16(IFGCPO(IGC)+1)
          ENDDO
        ELSE IF(NGREAC .GT. 0) THEN
*----
*  IF NGCOND = 0
*  AND NGREAC > 0
*  IFGCPO = IFGEDI
*----
          NGCCPO=NGREAC
          ENECPO(1)=ENET16(1)
          DO IGC=1,NGCCPO
            IFGCPO(IGC)=IFGEDI(IGC)
            ENECPO(IGC+1)=ENET16(IFGCPO(IGC)+1)
          ENDDO
        ELSE
*----
*  IF NGCOND = 0
*  AND NGREAC = 0
*  IFGCPO = IFGMTR
*----
          NGCCPO=NGMTR
          ENECPO(1)=ENET16(1)
          DO IGC=1,NGCCPO
            IFGCPO(IGC)=IFGMTR(IGC)
            ENECPO(IGC+1)=ENET16(IFGCPO(IGC)+1)   
          ENDDO
        ENDIF
      ENDIF
      IF(NGREAC .GT. 0) THEN
*----
*  IF NGREAC > 0
*    TEST IF CONDENSATION STRUCTURE PROVIDED BY IFGEDI
*    COMPATIBLE WITH IFGCPO AND IFGMTR
*  ENDIF
*----
        
        DO IGC=1,NGCCPO
          DO IGR=IGC,NGREAC
            IF(IFGCPO(IGC) .EQ. IFGEDI(IGR)) THEN
              IFGEDI(IGC)=IGR
              GO TO 135
            ENDIF
          ENDDO
          CALL XABORT(NAMSBR//
     >    ': IFGCPO AND IFGEDI NOT COMPATIBLE')
 135      CONTINUE
        ENDDO
      ENDIF
*----
*  NGMTR > 0
*    TEST IF CONDENSATION STRUCTURE PROVIDED BY IFGMTR
*    COMPATIBLE WITH IFGCPO
*  ENDIF
*----
      DO IGC=1,NGCCPO
        DO IGR=IGC,NGMTR
          IF(IFGCPO(IGC) .EQ. IFGMTR(IGR)) THEN
            IFGMTR(IGC)=IGR
            GO TO 155
          ENDIF
        ENDDO
        CALL XABORT(NAMSBR//
     >  ': IFGCPO AND IFGMTR NOT COMPATIBLE')
 155    CONTINUE
      ENDDO
      IF(IPRINT .GE.10) THEN
        IF(NGCOND .GT. 0)
     >    WRITE(IOUT,6022) (IFGCND(IGC),IGC=1,NGCOND)
        IF(NGCCPO .GT. 0) THEN
          WRITE(IOUT,6023) (IFGCPO(IGC),IGC=1,NGCCPO)
          WRITE(IOUT,6033) (ENECPO(IGC),IGC=1,NGCCPO+1)
        ENDIF
        WRITE(IOUT,6001)
      ENDIF
      DEALLOCATE(VELT16,VELEDI,IFGCPO)
      RETURN
*----
*  PRINT FORMAT
*----
 6000 FORMAT(1X,5('*'),' OUTPUT FROM ',A6,1X,5('*'))
 6001 FORMAT(1X,30('*'))
 6010 FORMAT(6X,'NUMBER OF LIBRARY GROUPS          = ',I10/
     >       6X,'NUMBER OF MAIN TRANSPORT GROUPS   = ',I10/
     >       6X,'NUMBER OF EDITING GROUPS          = ',I10/
     >       6X,'NUMBER OF CONDENSATION GROUPS     = ',I10/
     >       6X,'NUMBER OF CPO GROUPS              = ',I10)
 6020 FORMAT(6X,'MAIN TRANSPORT FEW GROUPS IDENTIFIER '/
     >10(2X,I6))
 6021 FORMAT(6X,'EDIT FEW GROUPS IDENTIFIER '/
     >10(2X,I6))
 6022 FORMAT(6X,'CONDENSATION FEW GROUPS IDENTIFIER '/
     >10(2X,I6))
 6023 FORMAT(6X,'CPO FEW GROUPS IDENTIFIER '/
     >10(2X,I6))
 6030 FORMAT(6X,'INITIAL ENERGY STRUCTURE (EV)'/
     >1P,10(2X,E10.3))
 6031 FORMAT(6X,'ENERGY STRUCTURE IN MAIN GROUPS (EV)'/
     >1P,10(2X,E10.3))
 6032 FORMAT(6X,'ENERGY STRUCTURE IN EDIT GROUPS (EV)'/
     >1P,10(2X,E10.3))
 6033 FORMAT(6X,'FINAL ENERGY STRUCTURE (EV)'/
     >1P,10(2X,E10.3))
 6040 FORMAT(6X,'AVERAGE VELOCITY IN INITIAL GROUPS (CM/S)'/
     >1P,10(2X,E10.3))
 6041 FORMAT(6X,'AVERAGE VELOCITY IN MAIN GROUPS (CM/S)'/
     >1P,10(2X,E10.3))
 6042 FORMAT(6X,'AVERAGE VELOCITY IN EDIT GROUPS (CM/S)'/
     >1P,10(2X,E10.3))
      END