summaryrefslogtreecommitdiff
path: root/Dragon/src/SPHMAC.f
blob: b53313e001c16aa00445ea5a736e706a4b116796 (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
*DECK SPHMAC
      SUBROUTINE SPHMAC(IPMACR,IPRINT,NMERGE,NALBP,NGCOND,ISCAT,NW,
     1 NIFISS,ILEAKS,VOLMER,FLXMER,SUNMER,SIGT,SIGW,DIFF,ZLEAK,OUTG,
     2 ALB,LFISS)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Recovery of the reference merged/condensed set of cross sections
* to be used by an SPH homogenization algorithm.
*
*Copyright:
* Copyright (C) 2011 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
* IPMACR  pointer to the condensed macrolib (L_MACROLIB signature).
* IPRINT  print flag (equal to 0 for no print).
* NMERGE  number of merged regions.
* NALBP   number of physical albedos.
* NGCOND  number of condensed groups.
* ISCAT   scattering anisotropy in the reference set of cross sections
*         (=1 isotropic in LAB; =2 linearly-anisotropic in LAB).
* NW      type of weighting for PN cross section info (=0 P0; =1 P1).
* NIFISS  number of fissile isotopes.
* ILEAKS  type of leakage calculation: =0 no leakage; =1 homogeneous
*         leakage (Diffon); =2 isotropic streaming (Ecco);
*         =3 anisotropic streaming (Tibere).
*
*Parameters: output
* VOLMER  merged volumes.
* FLXMER  merged/condensed averaged fluxes.
* SUNMER  merged/condensed production (fission + scattering) cross
*         sections. The third dimension is for secondary neutrons.
* SIGT    merged/condensed total P0 and P1 cross sections.
* SIGW    merged/condensed within-group scattering cross sections.
* DIFF    merged/condensed diffusion coefficients.
* ZLEAK   merged/condensed DB2 leakage rates.
* OUTG    merged/condensed leakage rates.
* ALB     physical albedos.
* LFISS   fission flag in mergeg zones.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPMACR
      INTEGER IPRINT,NMERGE,NALBP,NGCOND,ISCAT,NW,NIFISS,ILEAKS
      REAL VOLMER(NMERGE),FLXMER(NMERGE,NGCOND),
     1 SUNMER(NMERGE,NGCOND,NGCOND,ISCAT),SIGT(NMERGE,NGCOND,NW+1),
     2 SIGW(NMERGE,NGCOND,ISCAT+1),DIFF(NMERGE,NGCOND),
     3 ZLEAK(NMERGE,NGCOND),OUTG(NGCOND),ALB(NALBP,NGCOND)
      LOGICAL LFISS(NMERGE)
*----
*  LOCAL VARIABLES
*----
      PARAMETER(NSTATE=40)
      INTEGER ISTATE(NSTATE)
      DOUBLE PRECISION ZNUMER,ZLEAKA,ZDENUM
      CHARACTER HSIGN*12,SUFF*2,TEXT12*12
      TYPE(C_PTR) JPMACR,KPMACR
      INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,IPOS
      REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA,XSCAT,SIG1
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: PRODUC
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(NJJ(NMERGE),IJJ(NMERGE),IPOS(NMERGE))
      ALLOCATE(PRODUC(NMERGE,NGCOND,NIFISS),SIGMA(NMERGE*NIFISS),
     1 XSCAT(NMERGE*NGCOND))
*----
*  RECOVER MACROLIB INFORMATION
*----
      CALL LCMGTC(IPMACR,'SIGNATURE',12,HSIGN)
      IF(HSIGN.NE.'L_MACROLIB') CALL XABORT('SPHMAC: MACROLIB EXPECTED')
      CALL LCMGET(IPMACR,'STATE-VECTOR',ISTATE)
      IF(ISTATE(1).NE.NGCOND) CALL XABORT('SPHMAC: INVALID NGCOND')
      IF(ISTATE(2).NE.NMERGE) CALL XABORT('SPHMAC: INVALID NMERGE')
      NL=ISTATE(3)
      IF(ISTATE(4).NE.NIFISS) CALL XABORT('SPHMAC: INVALID NIFISS')
      IF(NIFISS.EQ.0) CALL XABORT('SPHMAC: NO FISSILE ZONES')
      ITRANC=ISTATE(6)
      ILEAK=ISTATE(9)
      IF(MAX(1,ISTATE(10)).NE.NW) CALL XABORT('SPHMAC: INVALID NW')
*----
*  SET OUTPUT INFORMATION TO ZERO
*----
      PRODUC(:NMERGE,:NGCOND,:NIFISS)=0.0
      FLXMER(:NMERGE,:NGCOND)=0.0
      SUNMER(:NMERGE,:NGCOND,:NGCOND,:ISCAT)=0.0
      SIGT(:NMERGE,:NGCOND,:NW+1)=0.0
      SIGW(:NMERGE,:NGCOND,:ISCAT+1)=0.0
      DIFF(:NMERGE,:NGCOND)=0.0
      ZLEAK(:NMERGE,:NGCOND)=0.0
      LFISS(:NMERGE)=.FALSE.
*----
*  RECOVER FLUX AND COMPUTE THE FISSION RATE INFORMATION
*----
      ZNUMER=0.0D0
      ZLEAKA=0.0D0
      ZDENUM=0.0D0
      CALL LCMGET(IPMACR,'VOLUME',VOLMER)
      IF(NALBP.GT.0) CALL LCMGET(IPMACR,'ALBEDO',ALB)
      JPMACR=LCMGID(IPMACR,'GROUP')
      DO 40 IGR=1,NGCOND
      KPMACR=LCMGIL(JPMACR,IGR)
      CALL LCMLEN(KPMACR,'FLUX-INTG',ILCMLN,ITYLCM)
      IF(ILCMLN.EQ.0) CALL XABORT('SPHMAC: MISSING FLUX-INTG INFO')
      CALL LCMGET(KPMACR,'FLUX-INTG',FLXMER(1,IGR))
      CALL LCMLEN(KPMACR,'NUSIGF',ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         CALL LCMGET(KPMACR,'NUSIGF',SIGMA(1))
         DO 35 IFIS=1,NIFISS
         DO 30 IBM=1,NMERGE
         SS=FLXMER(IBM,IGR)*SIGMA((IFIS-1)*NMERGE+IBM)
         PRODUC(IBM,IGR,IFIS)=PRODUC(IBM,IGR,IFIS)+SS
         ZNUMER=ZNUMER+SS
   30    CONTINUE
   35    CONTINUE
      ENDIF
   40 CONTINUE
*----
*  RECOVER EIGENVALUES
*----
      CALL LCMLEN(IPMACR,'K-EFFECTIVE',ILCMLN,ITYLCM)
      IF(ILCMLN.EQ.1) THEN
         CALL LCMGET(IPMACR,'K-EFFECTIVE',EIGENK)
      ELSE
         EIGENK=1.0
      ENDIF
      IF(IPRINT.GT.5) WRITE(6,'(/16H SPHMAC: EIGENK=,1P,E12.4)') EIGENK
      CALL LCMLEN(IPMACR,'B2  B1HOM',ILCMLN,ITYLCM)
      IF(ILCMLN.EQ.1) THEN
         CALL LCMGET(IPMACR,'B2  B1HOM',B2)
      ELSE
         B2=0.0
      ENDIF
      IF(IPRINT.GT.5) WRITE(6,'(/12H SPHMAC: B2=,1P,E12.4)') B2
*----
*  RECOVER MERGED/CONDENSED CROSS SECTIONS
*----
      DO 175 IGR=1,NGCOND
      KPMACR=LCMGIL(JPMACR,IGR)
      CALL LCMLEN(KPMACR,'NTOT0',ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         CALL LCMGET(KPMACR,'NTOT0',SIGMA(1))
      ELSE
         CALL XABORT('SPHMAC: MISSING NTOT0 INFO')
      ENDIF
      IF(ITRANC.NE.0) THEN
*        TRANSPORT CORRECTION.
         ALLOCATE(SIG1(NMERGE))
         CALL LCMGET(KPMACR,'TRANC',SIG1)
         DO 45 IBM=1,NMERGE
         SIGMA(IBM)=SIGMA(IBM)-SIG1(IBM)
   45    CONTINUE
         DEALLOCATE(SIG1)
      ENDIF
      DO 50 IBM=1,NMERGE
      ZDENUM=ZDENUM+SIGMA(IBM)*FLXMER(IBM,IGR)
   50 CONTINUE
      IF(ILEAKS.EQ.1) THEN
         CALL LCMLEN(KPMACR,'DIFF',ILCMLN,ITYLCM)
         IF(ILCMLN.EQ.0) CALL XABORT('SPHMAC: UNABLE TO RECOVER DIFF R'
     >   //'ECORDS IN THE MACROLIB OBJECT.')
         CALL LCMGET(KPMACR,'DIFF',DIFF(1,IGR))
      ENDIF
      DO 60 IBM=1,NMERGE
      IF(ILEAKS.EQ.1) THEN
         ZLEAK(IBM,IGR)=DIFF(IBM,IGR)*B2*FLXMER(IBM,IGR)
      ELSE IF(ILEAKS.GT.1) THEN
         CALL XABORT('SPHMAC: LEAKAGE MODEL NOT IMPLEMENTED')
      ELSE
         ZLEAK(IBM,IGR)=0.0
      ENDIF
      SIGT(IBM,IGR,1)=SIGMA(IBM)
      ZLEAKA=ZLEAKA+ZLEAK(IBM,IGR)
   60 CONTINUE
*
      CALL LCMLEN(KPMACR,'CHI',ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         DO 75 IFIS=1,NIFISS
         CALL LCMGET(KPMACR,'CHI',SIGMA(1))
         DO 70 IBM=1,NMERGE
         DO 65 JGR=1,NGCOND
         SS=SIGMA((IFIS-1)*NMERGE+IBM)*PRODUC(IBM,JGR,IFIS)/EIGENK
         IF(SS.NE.0.0) LFISS(IBM)=.TRUE.
         SUNMER(IBM,JGR,IGR,1)=SUNMER(IBM,JGR,IGR,1)+SS
   65    CONTINUE
   70    CONTINUE
   75    CONTINUE
      ENDIF
      DO 90 IW=2,MIN(NW+1,10)
      WRITE(TEXT12,'(4HNTOT,I1)') IW-1
      CALL LCMLEN(KPMACR,TEXT12,ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         CALL LCMGET(KPMACR,TEXT12,SIGMA(1))
         DO 80 IBM=1,NMERGE
         SIGT(IBM,IGR,IW)=SIGMA(IBM)
   80    CONTINUE
      ELSE
         DO 85 IBM=1,NMERGE
         SIGT(IBM,IGR,IW)=SIGT(IBM,IGR,1)
   85    CONTINUE
      ENDIF
   90 CONTINUE
*----
*  PROCESS SCATTERING INFORMATION
*----
      DO 170 INL=1,ISCAT
      WRITE(SUFF,'(I2.2)') INL-1
      CALL LCMLEN(KPMACR,'SIGW'//SUFF,ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         CALL LCMGET(KPMACR,'SIGW'//SUFF,SIGMA(1))
      ELSE
         SIGMA(:NMERGE)=0.0
      ENDIF
      IF((ITRANC.NE.0).AND.(INL.EQ.1)) THEN
*        TRANSPORT CORRECTION.
         ALLOCATE(SIG1(NMERGE))
         CALL LCMGET(KPMACR,'TRANC',SIG1)
         DO 120 IBM=1,NMERGE
         SIGMA(IBM)=SIGMA(IBM)-SIG1(IBM)
  120    CONTINUE
         DEALLOCATE(SIG1)
      ENDIF
      CALL LCMLEN(KPMACR,'NJJS'//SUFF,ILCMLN,ITYLCM)
      IF(ILCMLN.GT.0) THEN
         CALL LCMGET(KPMACR,'NJJS'//SUFF,NJJ)
         CALL LCMGET(KPMACR,'IJJS'//SUFF,IJJ)
         CALL LCMGET(KPMACR,'IPOS'//SUFF,IPOS)
         CALL LCMGET(KPMACR,'SCAT'//SUFF,XSCAT)
         DO 150 IBM=1,NMERGE
         IPO=IPOS(IBM)
         DO 130 JGR=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
         IF(IGR.EQ.JGR) THEN
            SS=SIGMA(IBM)*FLXMER(IBM,JGR)
            SIGW(IBM,IGR,INL)=SIGW(IBM,IGR,INL)+SS
         ELSE
            SS=XSCAT(IPO)*FLXMER(IBM,JGR)
         ENDIF
         SUNMER(IBM,JGR,IGR,INL)=SUNMER(IBM,JGR,IGR,INL)+SS
         IF(INL.EQ.1) ZDENUM=ZDENUM-SS
         IPO=IPO+1
  130    CONTINUE
  150    CONTINUE
      ENDIF
  170 CONTINUE
  175 CONTINUE
*----
*  COMPUTE REFERENCE LEAKAGE RATES
*----
      DO 182 IGR=1,NGCOND
      OUTG(IGR)=0.0
      DO 181 IBM=1,NMERGE
      OUTG(IGR)=OUTG(IGR)-SIGT(IBM,IGR,1)*FLXMER(IBM,IGR)-
     1 ZLEAK(IBM,IGR)
      DO 180 JGR=1,NGCOND
      OUTG(IGR)=OUTG(IGR)+SUNMER(IBM,JGR,IGR,1)
  180 CONTINUE
  181 CONTINUE
  182 CONTINUE
*
      DO 202 INL=1,ISCAT
      DO 201 IGR=1,NGCOND
      DO 200 IBM=1,NMERGE
      IF(VOLMER(IBM).NE.0.0) THEN
         SIGW(IBM,IGR,INL)=SIGW(IBM,IGR,INL)/FLXMER(IBM,IGR)
         DO 190 JGR=1,NGCOND
         SUNMER(IBM,JGR,IGR,INL)=SUNMER(IBM,JGR,IGR,INL)/FLXMER(IBM,JGR)
  190    CONTINUE
      ENDIF
  200 CONTINUE
  201 CONTINUE
  202 CONTINUE
      DO 215 IGR=1,NGCOND
      DO 210 IBM=1,NMERGE
      IF(VOLMER(IBM).NE.0.0) THEN
         ZLEAK(IBM,IGR)=ZLEAK(IBM,IGR)/FLXMER(IBM,IGR)
         FLXMER(IBM,IGR)=FLXMER(IBM,IGR)/VOLMER(IBM)
      ENDIF
  210 CONTINUE
  215 CONTINUE
*----
*  PRINT INFORMATION
*----
      IF(IPRINT.GT.4) THEN
         WRITE(6,'(/33H SPHMAC: type of PN weighting NW=,I2)') NW
         WRITE(6,240) ZNUMER/ZDENUM,ZNUMER/(ZDENUM+ZLEAKA)
         WRITE(6,250) 'VOLMER',(VOLMER(IKK),IKK=1,NMERGE)
         DO 220 IW=1,NW+1
         WRITE(TEXT12,'(4HNTOT,I1)') IW-1
         WRITE(6,250) TEXT12,((SIGT(IKK,IGR,IW),IKK=1,NMERGE),
     >   IGR=1,NGCOND)
  220    CONTINUE
         WRITE(6,250) 'FLXMER',((FLXMER(IKK,IGR),IKK=1,NMERGE),
     >   IGR=1,NGCOND)
         WRITE(6,250) 'ZLEAK',((ZLEAK(IKK,IGR),IKK=1,NMERGE),
     >   IGR=1,NGCOND)
         IF(NALBP.GT.0) THEN
           WRITE(6,250) 'ALBEDO',((ALB(IAL,IGR),IAL=1,NALBP),
     >     IGR=1,NGCOND)
           WRITE(6,250) 'OUTG',(OUTG(IGR),IGR=1,NGCOND)
         ENDIF
         DO 230 INL=1,ISCAT
         WRITE(SUFF,'(I2.2)') INL-1
         WRITE(6,250) 'SIGW'//SUFF,((SIGW(IKK,IGR,INL),IKK=1,NMERGE)
     >   ,IGR=1,NGCOND)
         WRITE(6,250) 'SUNMER'//SUFF,(((SUNMER(IKK,IGR,JGR,INL),
     >   IKK=1,NMERGE),IGR=1,NGCOND),JGR=1,NGCOND)
  230    CONTINUE
         IF(ILEAKS.EQ.1) THEN
            WRITE(6,250) 'DIFF',((DIFF(IKK,IGR),IKK=1,NMERGE),IGR=1,
     >      NGCOND)
         ENDIF
         WRITE(6,260) 'LFISS',(LFISS(IKK),IKK=1,NMERGE)
      ENDIF
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(XSCAT,SIGMA,PRODUC)
      DEALLOCATE(IPOS,IJJ,NJJ)
      RETURN
*
  240 FORMAT(/20H SPHMAC: K-INFINITY=,1P,D13.6/8X,12HK-EFFECTIVE=,D13.6,
     > 25H (FUNDAMENTAL MODE VALUE))
  250 FORMAT(/26H SPHMAC: VALUES OF VECTOR ,A,4H ARE/(1X,1P,10E13.5))
  260 FORMAT(/26H SPHMAC: VALUES OF VECTOR ,A,4H ARE/(1X,20L6))
      END