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
|
*DECK SHIEQU
SUBROUTINE SHIEQU(IPLIB,LEVEL,NGRO,NBISO,NBM,NBNRS,NRAT,MIX,
1 ISONAM,NOCONV,ISONR,GC,COEF,DENOM,XCOEF,XDENO,DEN,IMPX,SN)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the equivalent dilution.
*
*Copyright:
* Copyright (C) 2004 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
* IPLIB pointer to the internal microscopic cross section library
* (L_LIBRARY signature).
* LEVEL type of approximation (=1 use original Stamm'ler and
* Nordheim approximations; =2 use Riemann integration method
* with Nordheim approximation).
* NGRO number of energy groups.
* NBISO number of isotopes present in the calculation domain.
* NBM number of mixtures in the macrolib.
* NBNRS number of totally correlated resonant regions.
* NRAT number of terms in the rational pij expansion.
* MIX mix number of each isotope.
* ISONAM alias name of isotopes.
* NOCONV mixture convergence flag (.TRUE. if mixture IBM
* is not converged in group L).
* ISONR resonant isotope indices.
* GC Goldstein-Cohen parameters.
* COEF zone-independent weights.
* DENOM zone-independent base points.
* XCOEF zone-dependent weights.
* XDENO zone-dependent base points.
* DEN isotopic number density.
* IMPX print flag (equal to zero for no print).
*
*Parameters: output
* SN equivalent dilution.
*
*-----------------------------------------------------------------------
*
USE GANLIB
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPLIB
LOGICAL NOCONV(NBM,NGRO)
INTEGER LEVEL,NGRO,NBISO,NBM,NBNRS,NRAT,MIX(NBISO),
1 ISONAM(3,NBISO),ISONR(NBNRS),IMPX
REAL GC(NGRO,NBNRS),DEN(NBISO),SN(NGRO,NBISO)
COMPLEX COEF(NRAT,NGRO),DENOM(NRAT,NGRO)
COMPLEX*16 XCOEF(NRAT,NBNRS,NGRO),XDENO(NRAT,NBNRS,NGRO)
*----
* LOCAL VARIABLES
*----
PARAMETER(MAXRAT=9)
TYPE(C_PTR) JPLIB,KPLIB
LOGICAL LD
CHARACTER HNAMIS*12,HSMG*131
COMPLEX*16 EAV,TTT,CBS(MAXRAT)
INTEGER, ALLOCATABLE, DIMENSION(:) :: NFS
REAL, ALLOCATABLE, DIMENSION(:) :: EBIN,TBIN,SBIN
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(NFS(NGRO))
*
IF(NRAT.GT.MAXRAT) CALL XABORT('SHIEQU: MAXRAT OVERFLOW.')
JPLIB=LCMGID(IPLIB,'ISOTOPESLIST')
DO 200 IRS=1,NBNRS
ISO=ISONR(IRS)
KPLIB=LCMGIL(JPLIB,ISO) ! set ISO-th isotope
WRITE(HNAMIS,'(3A4)') (ISONAM(I0,ISO),I0=1,3)
CALL LCMLEN(KPLIB,'BIN-NFS',LENBIN,ITYLCM)
IF((LENBIN.GT.0).AND.(LEVEL.EQ.2)) THEN
CALL LCMGET(KPLIB,'BIN-NFS',NFS)
LBIN=0
DO 10 IGRP=1,NGRO
LBIN=LBIN+NFS(IGRP)
10 CONTINUE
ALLOCATE(EBIN(LBIN+1),TBIN(LBIN),SBIN(LBIN))
CALL LCMGET(KPLIB,'BIN-ENERGY',EBIN)
CALL LCMGET(KPLIB,'BIN-NTOT0',TBIN)
CALL LCMGET(KPLIB,'BIN-SIGS00',SBIN)
ELSE
NFS(:NGRO)=0
ENDIF
LBIN=0
DO 110 IGRP=1,NGRO
IF(NOCONV(MIX(ISO),IGRP)) THEN
EAV=0.0
IF(NFS(IGRP).GT.0) THEN
* RIEMANN INTEGRATION METHOD WITH NEWTON-RAPHSON ACCELERATION.
IF(IMPX.GE.10) THEN
WRITE(6,'(/17H SHIEQU: WEIGHTS:)')
WRITE(6,290) IGRP,(COEF(I,IGRP),I=1,NRAT)
WRITE(6,'(/21H SHIEQU: BASE POINTS:)')
WRITE(6,290) IGRP,(DENOM(I,IGRP),I=1,NRAT)
ENDIF
DO 20 I=1,NRAT
EAV=EAV+COEF(I,IGRP)*SQRT(DENOM(I,IGRP))
20 CONTINUE
SNI=(DBLE(EAV)**2)/DEN(ISO)
UG=0.0
BA=0.0
BS=0.0
DO 30 IGF=1,NFS(IGRP)
DELM=LOG(EBIN(LBIN+IGF)/EBIN(LBIN+IGF+1))
UG=UG+DELM
SIGA=MAX(0.0,TBIN(LBIN+IGF)-SBIN(LBIN+IGF))
SIGS=GC(IGRP,IRS)*MAX(0.0,SBIN(LBIN+IGF))
BA=BA+DELM*SIGA/(SIGA+SIGS+SNI)
BS=BS+DELM*SIGS/(SIGA+SIGS+SNI)
30 CONTINUE
BA=BA/UG
BS=BS/UG
ZCAL=SNI*BA/(1.0D0-BS)
DO 50 I=1,NRAT
CBS(I)=0.0
IF(XCOEF(I,IRS,IGRP).EQ.0.0) GO TO 50
EAV=XDENO(I,IRS,IGRP)/DEN(ISO)
DO 40 IGF=1,NFS(IGRP)
DELM=LOG(EBIN(LBIN+IGF)/EBIN(LBIN+IGF+1))
SIGA=MAX(0.0,TBIN(LBIN+IGF)-SBIN(LBIN+IGF))
SIGS=GC(IGRP,IRS)*MAX(0.0,SBIN(LBIN+IGF))
CBS(I)=CBS(I)+DELM*SIGS/(SIGA+SIGS+EAV)
40 CONTINUE
CBS(I)=CBS(I)/UG
CBS(I)=1.0D0/(1.0D0-CBS(I))
50 CONTINUE
TAUXA=0.0D0
DO 60 IGF=1,NFS(IGRP)
DELM=LOG(EBIN(LBIN+IGF)/EBIN(LBIN+IGF+1))
SIGA=MAX(0.0,TBIN(LBIN+IGF)-SBIN(LBIN+IGF))
SIGS=GC(IGRP,IRS)*MAX(0.0,SBIN(LBIN+IGF))
TTT=0.0D0
DO 55 I=1,NRAT
IF(XCOEF(I,IRS,IGRP).EQ.0.0) GO TO 55
EAV=XDENO(I,IRS,IGRP)/DEN(ISO)
TTT=TTT+XCOEF(I,IRS,IGRP)*EAV*CBS(I)/(SIGA+SIGS+EAV)
55 CONTINUE
TAUXA=TAUXA+DELM*SIGA*MAX(0.0D0,DBLE(TTT))
60 CONTINUE
IF(TAUXA.EQ.0.0) THEN
SNI=1.0E10
GO TO 90
ENDIF
TAUXA=TAUXA/UG
ITER=0
70 ITER=ITER+1
IF(IMPX.GE.10) THEN
WRITE(6,'(15H SHIEQU: GROUP=,I4,11H ITERATION=,I3,
1 10H DILUTION=,1P,E11.4,16H REFERENCE RATE=,E11.4,
2 13H APPROXIMATE=,E11.4)') IGRP,ITER,SNI,TAUXA,ZCAL
ENDIF
IF(ABS(TAUXA-ZCAL).LE.1.0D-5*ABS(TAUXA)) GO TO 90
IF(ITER.GT.20) THEN
WRITE(6,'(15H SHIEQU: GROUP=,I4,10H DILUTION=,1P,E11.4,
1 16H REFERENCE RATE=,E11.4,13H APPROXIMATE=,E11.4,
2 9H ISOTOPE=,A12,1H.)') IGRP,SNI,TAUXA,ZCAL,HNAMIS
IF(ABS(TAUXA-ZCAL).LE.5.0E-2*ABS(TAUXA)) THEN
WRITE(6,'(24H SHIEQU: *** WARNING ***)')
GO TO 90
ENDIF
CALL XABORT('SHIEQU: CONVERGENCE FAILURE.')
ENDIF
BA=0.0
BS=0.0
DBA=0.0
DBS=0.0
DO 80 IGF=1,NFS(IGRP)
DELM=LOG(EBIN(LBIN+IGF)/EBIN(LBIN+IGF+1))
SIGA=MAX(0.0,TBIN(LBIN+IGF)-SBIN(LBIN+IGF))
SIGS=GC(IGRP,IRS)*MAX(0.0,SBIN(LBIN+IGF))
BA=BA+DELM*SIGA/(SIGA+SIGS+SNI)
BS=BS+DELM*SIGS/(SIGA+SIGS+SNI)
DBA=DBA-DELM*SIGA/(SIGA+SIGS+SNI)**2
DBS=DBS-DELM*SIGS/(SIGA+SIGS+SNI)**2
80 CONTINUE
BA=BA/UG
BS=BS/UG
DBA=DBA/UG
DBS=DBS/UG
ZCAL=SNI*BA/(1.0D0-BS)
DZCAL=BA/(1.0D0-BS)+SNI*DBA/(1.0D0-BS)+SNI*BA*DBS/
1 (1.0D0-BS)**2
IF(DZCAL.LT.1.0D-15*ZCAL) GO TO 90
SNI=MAX(1.0D0,SNI-REAL((ZCAL-TAUXA)/DZCAL))
IF(SNI.GE.1.0E10) THEN
SNI=1.0E10
GO TO 90
ENDIF
GO TO 70
90 IF(IMPX.GE.5) THEN
WRITE(6,'(16H SHIEQU: REGION=,I3,7H GROUP=,I5,7H ISOTOP,
1 3HE='',A12,19H'' NB.OF ITERATIONS=,I3,1H.)') IRS,IGRP,
2 HNAMIS,ITER
ENDIF
ELSE
* ORIGINAL STAMM'LER APPROXIMATION.
DO 100 I=1,NRAT
EAV=EAV+XCOEF(I,IRS,IGRP)*SQRT(XDENO(I,IRS,IGRP))
100 CONTINUE
SNI=MAX(1.0D0,(DBLE(EAV)**2)/DEN(ISO))
ENDIF
SN(IGRP,ISO)=REAL(SNI)
IF(SN(IGRP,ISO).LE.0.0) THEN
WRITE (HSMG,300) HNAMIS,SN(IGRP,ISO),IGRP
CALL XABORT(HSMG)
ENDIF
ENDIF
LBIN=LBIN+NFS(IGRP)
110 CONTINUE
IF((LENBIN.GT.0).AND.(LEVEL.EQ.2)) THEN
DEALLOCATE(SBIN,TBIN,EBIN)
ENDIF
IF(IMPX.GE.5) THEN
LD=.FALSE.
DO 120 IGRP=1,NGRO
LD=LD.OR.NOCONV(MIX(ISO),IGRP)
120 CONTINUE
IF(LD) WRITE(6,310) HNAMIS,(SN(IGRP,ISO),IGRP=1,NGRO)
ENDIF
200 CONTINUE
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(NFS)
RETURN
*
290 FORMAT(10H GROUP NB.,I4,3X,1P,1H(,2E12.4,1H),:,2H (,2E12.4,1H),
1 2H (,2E12.4,1H),:,2H (,2E12.4,1H),:/(15X,1H(,2E12.4,1H),:,2H (,
2 2E12.4,1H),:,2H (,2E12.4,1H),:,2H (,2E12.4,1H)))
300 FORMAT(30HSHIEQU: THE RESONANT ISOTOPE ',A12,16H' HAS A NEGATIVE,
1 25H DILUTION CROSS-SECTION (,1P,E14.4,0P,10H) IN GROUP,I4,1H.)
310 FORMAT(/31H SHIEQU: DILUTIONS OF ISOTOPE ',A12,2H':/(1P,10E12.4))
END
|
