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
|
*DECK MOCFFIS
SUBROUTINE MOCFFIS(SUBSCH,NR,NS,NUN,MT,LINE,SEGLEN,NRSEG,NE,
1 MATALB,SIGANG,KEYFLX,YG,FLUX,EXPT,EXP2,FLM,
2 FLP,CYM,CYP,IDIR,OMG2)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solution of transport equation on a cyclic track
* ray-tracing (isotropic scattering case, 'source term isolation' off).
*
*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): R. Roy and R. Le Tellier
*
*Parameters: input
* SUBSCH track coefficients calculation subroutine.
* NR number of volumes.
* NS number of surfaces.
* NUN total number of unknowns in vectors FLUX.
* MT number of material mixtures.
* LINE number of segments on this tracking line.
* SEGLEN vector containing the lenght of the different segments of this
* track.
* NRSEG vector containing the region number of the different segments
* of this track.
* NE order of the polar quadrature set.
* MATALB index-number of the mixture assigned to each volume
* and the albedo to each surface.
* SIGANG total cross-sections and albedos.
* KEYFLX position of flux elements in FLUX vector.
* YG inverse of polar quadrature cosines.
* FLM total source vector for + direction.
* FLP total source vector for - direction.
* IDIR direction of fundamental current for TIBERE with MoC
* =0,1,2,3.
* OMG2 x, y and z components of the $3\\Omega^2$.
*
*Parameters: output
* FLUX vector containing the zonal flux moments.
*
*Parameters: scratch
* EXPT track coefficient.
* EXP2 quadratic expansion of (1-exp(-a*L))/L with small argument.
* CYM undefined.
* CYP undefined.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NR,NS,NUN,MT,LINE,NRSEG(LINE),NE,MATALB(-NS:NR),KEYFLX(NR)
INTEGER IDIR
REAL SIGANG(-6:MT),YG(NE)
DOUBLE PRECISION SEGLEN(LINE),FLUX(NUN),EXPT(NE,LINE),
1 EXP2(2,NE,LINE),FLM(NE,LINE),FLP(NE,LINE),CYM(NE,LINE),
2 CYP(NE,LINE),OMG2(NE,3)
EXTERNAL SUBSCH
*----
* LOCAL VARIABLES
*----
INTEGER MXE
PARAMETER (MXE=64)
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.D0,ZERO=0.D0)
INTEGER IE,NOIL,IL,JL,IND,INDC
DOUBLE PRECISION TEMP,TEMP2,FL0(MXE),CY0(MXE),FL1(MXE),CY1(MXE)
*
FL0(:NE)= ZERO
FL1(:NE)= ZERO
CY0(:NE)= ONE
CY1(:NE)= ONE
*----
* Calculation of the coefficients for this track.
*----
* MOCSCAS: Step-Characteristics Scheme with Tabulated Exponentials
* MOCDDFS: Diamond-Differencing Scheme
* MOCSCES: Step-Characteristics Scheme with Exact Exponentials
CALL SUBSCH(LINE,NR,NS,MT,NRSEG,MATALB,SEGLEN,SIGANG(-6),
1 EXPT,EXP2,NE,YG)
*----
* Summation along the track in both directions
*----
DO 4226 IL= 1, LINE
JL= LINE + 1 - IL
DO 4225 IE=1,NE
* phi_k
TEMP = FL0(IE)
TEMP2 = FLM(IE,IL)
* phi_{k+1}
FL0(IE) = TEMP*EXPT(IE,IL) + TEMP2*EXP2(1,IE,IL)
* phi_k * ((1 - exp(-tau_k)) / tau_k)
FLM(IE,IL) = TEMP*EXP2(1,IE,IL) + TEMP2*EXP2(2,IE,IL)
* ((1 - exp(-tau_k)) / tau_k) * exp(-tau_1^{k-1})
CYM(IE,IL) = CY0(IE) * EXP2(1,IE,IL)
* exp(-tau_1^{k})
CY0(IE) = CY0(IE) * EXPT(IE,IL)
*
TEMP = FL1(IE)
TEMP2 = FLP(IE,JL)
FL1(IE) = TEMP*EXPT(IE,JL) + TEMP2*EXP2(1,IE,JL)
FLP(IE,JL) = TEMP*EXP2(1,IE,JL) + TEMP2*EXP2(2,IE,JL)
CYP(IE,JL) = CY1(IE) * EXP2(1,IE,JL)
CY1(IE) = CY1(IE) * EXPT(IE,JL)
4225 CONTINUE
4226 CONTINUE
DO 4230 IE=1,NE
TEMP=ONE-CY0(IE)
IF (TEMP.GT.ZERO) THEN
FL0(IE)= FL0(IE)/TEMP
ELSE
FL0(IE)= ZERO
ENDIF
TEMP=ONE-CY1(IE)
IF (TEMP.GT.ZERO) THEN
FL1(IE)= FL1(IE)/TEMP
ELSE
FL1(IE)= ZERO
ENDIF
4230 CONTINUE
DO 4240 IL= 1, LINE
NOIL = NRSEG(IL)
IF( NOIL.GT.0 )THEN
IND=KEYFLX(NOIL)
INDC=NUN/2+IND
DO 4241 IE=1,NE
FLUX(IND)= FLUX(IND)
> + ((FL0(IE)*CYM(IE,IL)+FLM(IE,IL))
> +(FL1(IE)*CYP(IE,IL)+FLP(IE,IL)))
* CALCULATE XI, YI OR ZI FOR TIBERE
IF(IDIR.GT.0) THEN
FLUX(INDC)= FLUX(INDC)
> + ((FL0(IE)*CYM(IE,IL)+FLM(IE,IL))
> +(FL1(IE)*CYP(IE,IL)+FLP(IE,IL)))
> *OMG2(IE,IDIR)
ENDIF
4241 CONTINUE
ENDIF
4240 CONTINUE
*
RETURN
END
|
