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
|
*DECK LIBCOM
SUBROUTINE LIBCOM(NFS,DELTA,SIGAF,SIGTF,NORA,NOR,COMOM)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute a set of comoments ((SIGA**P)*(SIGT**Q)).
*
*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): A. Hebert
*
*Parameters: input
* NFS number of fine energy groups.
* DELTA lethargy widths of the fine groups.
* SIGAF microscopic absorption x-sections in the fine groups.
* SIGTF microscopic total x-sections in the fine groups.
* NORA related to the number of absorption moments to preserve.
* NOR related to the number of total moments to preserve:
* (2-NORA)/2 <= P <= (NORA+1)/2 and (2-NOR)/2 <= Q <= (NOR+1)/2.
*
*Parameters: output
* COMOM comoments.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NFS,NORA,NOR
REAL DELTA(NFS),SIGAF(NFS),SIGTF(NFS)
DOUBLE PRECISION COMOM(NORA,NOR)
*----
* LOCAL VARIABLES
*----
DOUBLE PRECISION DEL,T,T0,SIGT,SIGA
INTEGER PNOR,QNOR
*
DEL=0.0D0
DO 15 PNOR=1,NORA
DO 10 QNOR=1,NOR
COMOM(PNOR,QNOR)=0.0D0
10 CONTINUE
15 CONTINUE
*
DO 80 IGF=1,NFS
DELF=DELTA(IGF)
SIGT=MAX(0.001,SIGTF(IGF))
SIGA=MAX(0.001,SIGAF(IGF))
DEL=DEL+DELF
T0=DELF
DO 40 PNOR=MAX(1,NORA/2),NORA
T=T0
DO 20 QNOR=MAX(1,NOR/2),NOR
COMOM(PNOR,QNOR)=COMOM(PNOR,QNOR)+T
T=T*SIGT
20 CONTINUE
T=T0/SIGT
DO 30 QNOR=NOR/2-1,1,-1
COMOM(PNOR,QNOR)=COMOM(PNOR,QNOR)+T
T=T/SIGT
30 CONTINUE
T0=T0*SIGA
40 CONTINUE
T0=DELF/SIGA
DO 70 PNOR=NORA/2-1,1,-1
T=T0
DO 50 QNOR=MAX(1,NOR/2),NOR
COMOM(PNOR,QNOR)=COMOM(PNOR,QNOR)+T
T=T*SIGT
50 CONTINUE
T=T0/SIGT
DO 60 QNOR=NOR/2-1,1,-1
COMOM(PNOR,QNOR)=COMOM(PNOR,QNOR)+T
T=T/SIGT
60 CONTINUE
T0=T0/SIGA
70 CONTINUE
80 CONTINUE
*
IF(DEL.EQ.0.0) CALL XABORT('LIBCOM: ALGORITHM FAILURE.')
DO 100 PNOR=1,NORA
DO 90 QNOR=1,NOR
COMOM(PNOR,QNOR)=COMOM(PNOR,QNOR)/DEL
90 CONTINUE
100 CONTINUE
RETURN
END
|
