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
|
*DECK AFMTAV
SUBROUTINE AFMTAV (NBURN,ITM,XBMAX,XBMIN,YS,NBMIN,NBMAX,XB,SIGAV)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Time average calculation using different approximation.
*
*Copyright:
* Copyright (C) 1996 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):
* M.T. Sissaoui
*
*Parameters: input
* NBURN total number of steps.
* ITM type of the approximation (1-Lagrange; 2-spline; 3-Hermite)
* XBMAX highest value.
* XBMIN lower value.
* YS parameter to be integrated
* NBMIN
* NBMAX
* XB steps
*
*Parameters: output
* SIGAV average value of YS
*
*-----------------------------------------------------------------------
*
REAL YS(NBURN),XB(NBURN),SIGAV
REAL UU(2)
DOUBLE PRECISION DD
REAL, ALLOCATABLE, DIMENSION(:) :: Y,U
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(Y(NBURN),U(NBURN))
*
IF(NBMAX.GT.0) THEN
INMAX=NBMAX-2
INMIN=NBMIN
ELSE
INMAX=ABS(NBMAX)-1
INMIN=NBMIN
IF(ITM.EQ.1)
1 CALL XABORT('AFMTAV: MORE BURNUP STEPS ARE REQUIRED TO USE '
1 //' LAGRANGE METHOD, CHOOSE HERMIT OR SPLINE METHOD')
ENDIF
IF(ABS(NBMAX).GT.NBMIN) THEN
SIGAV=0.0
*
IF(ITM.EQ.1) THEN
* TIME AVERAGE CALCULATION USING LAGRANGE APPROXIMATION.
*
DO 113 IR=INMIN,INMAX
I1=IR
I2=IR+1
I3=IR+2
XBI=MAX(XBMIN,XB(IR))
XBF=MIN(XBMAX,XB(IR+1))
TX=XBF-XBI
TX2=XBF**2-XBI**2
TX3=XBF**3-XBI**3
X12=XB(I1)-XB(I2)
X13=XB(I1)-XB(I3)
X23=XB(I2)-XB(I3)
XA12=XB(I1)+XB(I2)
XA13=XB(I1)+XB(I3)
XA23=XB(I2)+XB(I3)
XM12=XB(I1)*XB(I2)
XM13=XB(I1)*XB(I3)
XM23=XB(I2)*XB(I3)
Y1=YS(I1)/(X12*X13)
Y2=-YS(I2)/(X12*X23)
Y3=YS(I3)/(X13*X23)
*
SIGAV=SIGAV +
1 Y1*(TX3/3.0-XA23*TX2/2.0+XM23*TX)+
1 Y2*(TX3/3.0-XA13*TX2/2.0+XM13*TX)+
1 Y3*(TX3/3.0-XA12*TX2/2.0+XM12*TX)
113 CONTINUE
*
ELSE IF(ITM.EQ.2) THEN
* TIME AVERAGE CALCULATION USING SPLINE APPROXIMATION.
* THE LOWER BOUNDARY CONDITION IS SET TO BE NATURAL
Y(1)=0.0
U(1)=0.0
* THE UPPER BOUNDARY CONDITION IS SET EITHER TO BE NATURAL
QN=0.0
UN=0.0
*
DO 103 IR=2,NBURN-1
SIG=(XB(IR)-XB(IR-1))/(XB(IR+1)-XB(IR-1))
P=SIG*Y(IR-1)+2.0
Y(IR)=(SIG-1.0)/P
U(IR)=(6.*((YS(IR+1)-YS(IR))/(XB(IR+1)-XB(IR))-
1 (YS(IR)-YS(IR-1))/(XB(IR)-XB(IR-1)))/(XB(IR+1)-
1 XB(IR-1))-SIG*U(IR-1))/P
103 CONTINUE
*
Y(NBURN)=(UN-QN*U(NBURN-1))/(QN*Y(NBURN-1)+1.0)
*
DO 104 K=NBURN-1,1,-1
Y(K)=Y(K)*Y(K+1)+U(K)
104 CONTINUE
*
* COMPUTE THE INTEGRAL OF THE X-SECTION
INMAX=NBMAX-2
INMIN=NBMIN
DO 300 IR=INMIN,INMAX
H=XB(IR+1)-XB(IR)
XBI=MAX(XBMIN,XB(IR))
XBF=MIN(XBMAX,XB(IR+1))
*
DB=XBF-XBI
HF=XB(IR+1)-XBF
HI=XB(IR+1)-XBI
*
AI=-0.5*(HF**2-HI**2)/H
BI=DB-AI
CI=-(AI/6)*H**2-(HF**4-HI**4)/(24*H)
DI=-(BI/6)*H**2-(HF**4-HI**4)/(24*H)
*
SIGAV=SIGAV+AI*YS(IR)+BI*YS(IR+1)+
1 CI*Y(IR)+DI*Y(IR+1)
300 CONTINUE
ELSE IF(ITM.EQ.3) THEN
* TIME AVERAGE CALCULATION USING HERMIT APPROXIMATION.
DO 101 I=1,NBURN
Y(I)=YS(I)
101 CONTINUE
* TAKE THE DERIVATIVE WITH RESPECT TO BURNUP OR NEUTRON EXPOSURE AT
* TABULATION POINTS.
CALL ALDERV(NBURN,XB,Y)
*
* COMPUTE THE INTEGRAL OF THE X-SECTION
DD=0.0D0
DO 200 IR=1,NBURN-1
IF((XBMIN.LT.XB(IR+1)).AND.(XBMAX.GT.XB(IR))) THEN
DX=XB(IR+1)-XB(IR)
XBI=MAX(XBMIN,XB(IR))
XBF=MIN(XBMAX,XB(IR+1))
CC=0.5*(XBF-XBI)
U1=(XBI-0.5*(XB(IR)+XB(IR+1)))/DX
U2=(XBF-0.5*(XB(IR)+XB(IR+1)))/DX
UU(1)=0.5*(-(U2-U1)*0.577350269189626+U1+U2)
UU(2)=0.5*((U2-U1)*0.577350269189626+U1+U2)
DO 190 J=1,2
H1=3.0*(0.5-UU(J))**2-2.0*(0.5-UU(J))**3
H2=(0.5-UU(J))**2-(0.5-UU(J))**3
H3=3.0*(0.5+UU(J))**2-2.0*(0.5+UU(J))**3
H4=-(0.5+UU(J))**2+(0.5+UU(J))**3
DD=DD+(H1*YS(IR)+H2*Y(IR)*DX+H3*YS(IR+1)+
1 H4*Y(IR+1)*DX)*CC
190 CONTINUE
ENDIF
200 CONTINUE
SIGAV=REAL(DD)
ENDIF
SIGAV=SIGAV/(XBMAX-XBMIN)
ELSE
SIGAV=YS(NBMIN)
ENDIF
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(U,Y)
RETURN
END
|
